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Preface 

This  research  began  on  two  fronts.  One  was  a  desire  to  build  upon  the  work 
performed  by  Lt  Carl  Tong  in  the  area  of  multisensor  data  fusion.  Lt  Tong  left  a 
variety  of  software  tools  which  proved  useful  primarily  in  displaying  images  on  the 
Evans  and  Sv- '  '.and  PS340  Graphics  Workstation.  The  second  was  an  interest  in 
creating  and  implementing  computer-generated  interferograms.  The  topic  of 
Hough  transforms  was  originally  chosen  as  an  application  for  computer-generated 
interferograms.  When  difficulties  in  using  the  interferograms  on  complex  fields 
arose,  discrete  implementations  and  properties  of  the  Hough  transform  became  the 
primary  focus. 

This  project  created  some  useful  software  tools  for  implementing  and  per¬ 
forming  digital  filtering  on  Hough  transforms  for  use  on  the  VAX  computer.  It 
also  implemented  an  ”in  house"  technique  for  producing  computer-generated  in¬ 
terferograms  which  proved  invaluable  to  several  other  research  projects. 

There  are  a  few  individuals  whom  I  would  like  to  acknowledge  for  their  assis¬ 
tance  throughout  my  thesis  effort.  First  and  foremost  is  my  thesis  advisor  Capt 
Steven  K.  Rogers  whose  continual  patience,  enthusiasm  and  prodding  gave  me  the 
necessary  direction  and  focus  for  my  research.  I  would  also  like  to  thank  my 
committee  members  Dr.  Mathew  Kabrisky  and  LTC  James  Mills  for  providing 
intellectual  stimulation.  Finally,  I  would  like  to  thank  fellow  AFTT  student  Lt  Mike 
Mayo  for  his  photographic  assistance. 


—  Marvin  L.  Hill 
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Abstract 

\j  This  thesis  applied  the  normal  straight  line  parameterization  of  the  Hough 
transform  to  a  variety  of  images  using  the  accumulator  method.  Simple  inputs 
were  used  initially  to  illustrate  the  distortion  characteristics  of  the  Hough  transform 
due  to  rotations,  scales  and  translations  of  an  input.  Making  use  of  work  per¬ 
formed  by  D.  Casasent  and  R.  Krishnapuram  of  Camegie-Mellon  University,  the 
Hough  transform  was  then  applied  to  segmented  and  edged  doppler  images.  A 
distort-and-compare  routine,  which  makes  use  of  the  Hough  space  distortion  char¬ 
acteristics,  was  implemented  in  the  Hough  transform  domain  to  estimate  input 
space  characteristics  -■'f  an  object. 

Next,  the  Hough  transform,  generated  using  Fourier  transform  techniques, 
was  applied  to  some  of  the  same  inputs  to  demonstrate  that  the  accumulator 
method  is  actually  a  discrete  version  of  the  continuous  Hough  transform.  An 
unsuccessful  attempt  at  implementing  the  continuous  Hough  transform  using  com¬ 
puter-generated  interferograms  was  outlined.  A  method  of  implementing  the  con¬ 
tinuous  Hough  transform  and  its  inverse  using  phase  filters  was  presented  as  a 
suggestion  for  further  research.  '  .  '  ^  ^ 


1.  Introduction 


A.  Background.  Mission  work  load  in  all  combat  areas  has  increased 
dramatically  in  recent  years.  A  significant  portion  of  combat  requirements 
involves  acquisition  and  identification  of  hostile  targets  within  the  threat 
environment.  Thus,  a  considerable  amount  of  research  has  been  devoted  to 
segmenting  targets  from  a  variety  of  backgrounds  and  current  algorithms  are  often 
quite  effective  in  performing  this  task.  The  ability  to  classify  targets  either  before 
or  during  segmentation  is  a  much  more  difficult  task.  Investigation,  analysis  and 
implementation  of  alternative  feature  extraction  techniques  is  essential  to  ensure 
out-numbered  forces  maintain  a  technological  edge.  The  Hough  transform  has 
shown  some  promise  in  its  feature  extraction  properties.  This  thesis  will 
investigate  the  properties  of  the  Hough  transform  and  research  methods  for  its 
implementation. 

B.  Approach  and  Scope.  The  first  task  in  this  thesis  project  is  to  gain  good 
working  knowledge  of  the  ADA  programming  language  and  programming 
techniques  in  general.  A  variety  of  display  drivers  and  input/output  routines  were 
created  in  the  Tong  [20]  thesis  project  with  most  of  the  programming 
accomplished  in  the  ADA  language.  Thus,  creating  the  necessary  requirement  of 
understanding  current  file  formatting  in  order  to  obtain,  create  or  process  image 
files.  The  second  task  is  to  research  current  knowledge  of  Hough  transform 
techniques  which  might  suggest  efficient  means  of  implementation.  Third, 
working  edge  images  must  be  created  to  test  Hough  transform  implementations 
and  determining  Hough  transform  characteristics.  Elementary  test  patterns  are 
required  as  well  as  examples  of  actual  segmented  data.  In  addition,  computer 
implementations  of  the  Hough  transform  techniques  must  be  created  and  display 


routines  modified  to  present  the  output.  The  final  task  will  include  an  investigation 
into  possible  real-time  implementations  of  the  Hough  transform. 

C.  Organization.  This  thesis  begins  with  a  summary  of  current  theory  in  Section 
n.  A  basic  definition  of  the  Hough  transform  is  presented  followed  by  some  of  the 
methods  used  for  its  implementation.  Current  attempts  to  obtain  the  transform 
optically  are  also  discussed.  The  creation  of  working  edge  images  is  briefly 
outlined  in  Section  in  followed  by  a  description  of  the  two  Hough  transform 
processes  in  Section  IV.  The  effects  of  the  transform  on  basic  shapes  is  analyzed 
in  this  section  as  well.  A  process  for  estimating  input  space  distortions  within  the 
Hough  domain  is  described  in  Section  V  and  examples  of  its  application  to  real 
edge  images  are  depicted.  This  is  followed  by  a  general  discussion  of  the  results 
of  the  thesis  in  Section  VI.  In  Section  Vn,  a  suggested  continuous  optical 
implementation  technique  is  presented  as  a  recommendation  for  further  research 
with  background  data  for  this  approach  outlined  in  Appendix  B  and  Appendix  C. 
Finally,  concluding  statements  on  the  research  project  reside  in  Section  VDI. 


II.  Sununaty  of  Current  Knowledge 


The  Hough  transform  was  developed  in  1962  by  P.V.C.  Hough  [1]  to  trace 
the  paths  of  atomic  particles  in  cloud  chambers.  The  initial  procedures  outlined  by 
Hough  [1]  were  written  in  terms  of  electrical  circuit  diagrams.  It  wasn’t  until  1972 
(Duda-Hart  [2])  that  Hough  transform  techniques  were  applied  to  pattern 
recognition. 

A.  Basic  Approach.  Current  literature  applies  a  rather  broad  definition  to  the 
term  Hough  transform.  In  general,  a  Hough  transform  maps  points  in  an  input 
feature  space  (e.g.,  a  Cartesian  coordinate  space)  to  curve  shape  parameters  in  an 
output  space.  This  mapping  is  usually  based  on  the  coordinates  of  the  point  and 
possibly  the  directional  gradient  at  that  point. 

Hough  transform  techniques  are  used  to  detect  specific  types  of  curves  within 
the  input  space.  Thus,  the  curve  type  is  normally  included  when  describing  a 
specific  type  of  Hough  transform.  Variations  most  often  discussed  include  Hough 
transform  mappings  to  straight  line,  circle  and  ellipsoid  parameter  spaces.  These 
variations,  along  with  generalized  parameter  space  applications,  are  discussed  in 
the  following  paragraphs  of  this  section. 

Implementation  of  Hough  transforms  for  detecting  curves  within  an  input 
scene  requires  segmentation  and  edging  of  that  input  scene.  This  is  not  an  overly 
restrictive  condition  since  many  segmentation  algorithms  have  been  developed 
(e.g.,  Tong  [20]  and  Ruck  [23]).  The  Hough  transform  can  then  be  applied  to 
e.xtract  features  from  the  outline  of  the  segmented  objects  as  an  aid  to 
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B.  Simple  Curves.  The  applications  of  the  Hough  transform  presented  by 
Duda-Hart  [2],  [3]  and  Shapiro  [4]  describe  mappings  to  a  straight  line  parameter 
space.  From  the  basic  slope  intercept  equation  for  a  line  in  Cartesian  coordinates, 

y  =  ax  +  b  (2.1) 

a  binary  input  space  f(x,y)  is  mapped  to  a  parameter  space  H(a,b).  Thus  all  points 
in  the  input  space  situated  along  the  line  y  =  aox  +  bo  are  mapped  to  a  single  point 
H(ao,bo)  in  the  Hough  space. 

Similarly,  a  Hough  transform  parameter  space  developed  for  detection  of 
paraboloids  of  the  form 

y*cx2  +  dx  +  e  (2.2) 

would  produce  a  Hough  output  space  H(c,d,e).  The  methods  used  to  implement 
this  type  of  Hough  transform  vary  considerably.  Most,  however,  use  a  variation  of 
the  procedure  used  by  Shapiro  [5:129]  which  require  taking  derivatives  at  each 
edge  point  (pixel)  in  the  input  space. 

C.  Hough  Transform  Using  a  R-Table.  The  first  useful  generalized  Hough 
transform  technique  was  introduced  by  Ballard  [6]  in  1981.  In  addition  to 
mappings  to  parameter  spaces  of  analytic  curves,  this  procedure  can  create  a 


Hough  parameter  space  for  an  arbitrary  shape.  The  transform  is  applied  by  taking 
directional  derivatives  at  each  edge  point.  This  value  is  compared  to  a  table 
describing  which  points  in  the  output  space  to  increment.  Each  point  in  the  output 
Hough  space  is  then  incremented  as  shown  in  Figure  2.1, 


Figure  2.1.  Geometry  for  Generalized  Hough  Transform[6:116] 

A  directional  gradient  is  taken  at  each  edge  point  x  in 
the  input  space.  The  gradient  value  corresponds  to  an 
increment  vector  r  in  the  R-Table.  The  increment  vector 
identifies  the  point,  a,  to  increment  in  the  Hough  space. 

with  the  edge  point  in  the  input  space,  directional  gradient,  increment  vector  and 
increment  point  in  the  output  space  are  represented  by  x,  <j>,  r  and  a  respectively. 
The  table,  referred  to  as  an  R-Table,  contains  a  set  of  vectors  indexed  to  each 
directional  derivative  of  the  curve  to  be  detected.  Rotations  and  scales  of  the  cur\e 
in  the  input  plane  can  be  handled  by  rotating  and  scaling  the  increment  vectors 
within  the  R-Table. 


D.  Normal  Straight  Line  Parameterization.  The  Hough  space  parameterization 
based  on  the  normal  representation  of  a  straight  line. 

p  »  X  cos(0)  +  y  sin(6)  (2.3) 

was  used  initially  to  alleviate  the  problems  associated  with  infinite  slopes  in  the 
alternative  slope-intercept  parameterization  of  Eqn  (2.1).  The  geometry  of  this 
parameterization  is  shown  in  Figure  2.2. 


One  of  the  major  advantages  of  this  Hough  space  parametenzation  is 
computational  simplicity.  Using  an  accumulator  technique,  the  normal  straight 
line  Hough  transform  is  quickly  computed  by  incrementing  points  in  the  Hough 
parameter  space  along  the  sinusoid  defined  by  Eqn  (2.3)  for  each  edge  point  (x,y) 
within  the  input  scene.  When  edge  points  in  an  input  scene  are  situated  along  a 
straight  line,  the  corresponding  sinusoids  will  intersect  at  a  single  point  causing  a 
ma.xima  at  that  point  within  the  Hough  space  (Figure  2.3). 


Casasent-Krishnapuram  [8],  [9]  documented  some  extremely  significant 

characteristics  of  the  normal  straight  line  Hough  space  parameterization. 

In  these  articles,  Casasent  and  Krishnapuram  studied  shifts,  rotations  and 
scales  of  edge  curves  within  input  scenes  and  showed  that  the  corresponding 
distortions  within  the  Hough  parameter  space  followed  well  defined  rules.  The 
relationship  between  the  Hough  transform  of  an  input  curve  and  the  Hough 
transform  of  its  scaled  version  was  shown  to  be 

Hs(  0,  p  )  =  H(  0.  s  p  )  (2.4) 

Thus,  scaling  of  the  input  space  by  s,  causes  the  Hough  space  to  be  scaled  along 
the  p  axis  [9:304]  by  the  same  value. 

An  object  rotated  by  4)  in  the  input  space  caused  the  relationship 

Hr(0.p)«H(0-4».p)  (2.5) 

where  the  original  Hough  space  is  shifted  along  0  by  4)  [9:305]. 

Translations  of  an  object  in  the  input  scene  such  that 

ft(  X,  y  )  =  f(  X  -  a,  y  -  b  ) 

creates  a  sinusoidal  distortion 


H, 


H(  0,  p  -  t  cos(  0  -  o( ), 


p  +  t  cos(0  -  a  )  >  0 


(2.6) 


where  t  =  ^  +  h?  and  a  = 

along  the  p  axis  of  the  Hough  space  [9:306]. 

The  two  cases  in  Eqn  (2.6)  are  merely  the  result  of  discarding  negative  values 
of  p  and  does  not  point  to  any  inherent  discontinuity  in  the  normal  straight  line 
parameterization.  If  the  redundant  representation  of  the  Hough  space  where  used, 
the  first  case  of  Eqn  (2.6)  would  be  sufficient  to  describe  the  Hough  space 
distortion. 

Combining  Eqns  (2.4),  (2.5)  and  (2.6)  gives  the  generalized  characterization 
for  distortion  of  the  Hough  parameter  space  for  scales,  rotations  and  shifts 
(translations)  of  an  object  in  the  input  space  [9:306]: 


f  H(  0  -  4>,  s  {  p  +  t  cos(0  -  a)  }  ),  P  +  t  cos(0  -  oi(  )  ^  0 

H’=  {  (2.7) 

I  H(0-<t)-Tr,s{-p-t  cos(0  -  a) } ),  p  +  t  cos(0  -  a  )  <  0 

Thus,  a  search  for  any  two-dimensional  object  within  an  input  scene  can  be 
performed  in  the  Hough  parameter  space  by  creating  a  Hough  transform  template 
and  using  a  distort  and  compare  routine  to  find  it’s  location,  rotation  and/or  scale. 


E.  Optical  Implementations.  The  normal  straight  line  Hough  transform  can  be 
shown  to  be  a  special  case  of  the  Radon  transform  (Appendix  A),  where 


H(  0,  p )  =  jy  f(  X,  y  )  8[  p  -  xcos(6)  -  y  sin(9)  ]  dx  dy  (2.8' 

-OO 

and 

f(  X,  y  )  =  1,  (  X,  y  )  is  an  edge  pixel 

0,  (  X,  y  )  is  not  an  edge  pixel 

Deans  [14:96-100]  describes  a  direct  relationship  between  the  Radon  transform 
and  the  Fourier  transform.  This  relationship  suggests  the  Radon  transform,  and 
thus  the  normal  straight  line  Hough  transform,  can  be  obtained  in  real  time  using 
some  type  of  optical  implementation. 

Eichmann  and  Dong  [13]  developed  a  method  of  optically  generating  slices 
of  the  Hough  transform.  As  shown  in  Figure  2.4,  the  input  image  is  placed  in 
plane  P-\,  a  slit  placed  in  plane  P2  passes  an  angular  slice  of  the  Fourier  transform 
of  the  input.  Plane  /’3  records  H(0i,  p)  for  each  rotation,  0j,  of  the  image  in  plane 
Pi. 

Other  optical  implementations  follow  the  same  form  as  the  Eichmann-Dong 
method.  Steier  and  Short  [11]  employed  a  similar  method.  They  used  a  rotating 
dove  prism  to  effectively  rotate  an  input  wavefront  and  a  cylindrical  lens  to 
eliminate  the  slit  in  plane  P2  of  Figure  2.4.  Another,  somewhat  more  elaborate, 
approach  was  used  by  Ambs  et  al  [12].  They  were  able  to  obtain  several  samples 
of  the  the  Hough  transform  through  a  Fourier  transform  filter  composed  of  a 
matrix  of  holograms.  No  method,  however,  has  been  developed  to  generate  a 
continuous  optical  reproduction  of  the  Hough  transform. 
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III.  Creation  of  Edge  Images 


A.  Source  Data.  The  edge  images  used  in  this  thesis  were  generated  from 
doppler  and  laser  radar  digital  image  data  obtained  from  the  Air  Force  Wright 
Aeronautical  Laboratories  (AFWAL)  Avionics  Laboratory.  Computational 
processing  of  this  data  was  performed  on  a  VAX  11/780  and  an  Evans  and 
Sutherland  PS340  Graphics  Workstation  was  used  for  displaying  images. 


B.  Edging  Process.  The  block  diagram  in  Figure  3.1  outlines  process  used  to 
create  edge  images  from  the  raw  doppler  and  laser  radar  data. 


Figure  3.1.  Block  Diagram  of  Edging  Process 


Raw  image  data  was  passed  through  a  Tong  [20:24-37]  gradient  magnitude  filter. 
The  result  was  thresholded  and  binarized.  Several  stages  of  median  filtering 
[20:49]  were  applied  to  remove  unwanted  noise.  This  results  in  a  solid  mask  of 
the  object.  Mask  is  then  sent  through  another  gradient  magnitude  filter  to  produce 
me  edge  images.  A  raw  doppler  image  and  its  corresponding  edge  image  are 
shown  in  Figures  3.2  and  3.3.  In  general,  the  full  Tong  algorithm  [21]  or  a  similar 
procedure  would  be  used  to  segment  the  mask  of  the  object  from  an  input  scene. 


Once  a  set  of  edge  images  has  been  created.  Hough  transform  processing  can 
begin.  The  next  chapter  discusses  two  methods  of  producing  the  Hough  transform 
along  with  some  examples  of  transformed  images  to  illustrate  the  effects  of 
distortions  of  an  input  object  on  the  Hough  transform  domain. 


rv.  Hough  Transform  Images 


Throughout  the  rest  of  this  thesis,  the  term  Hough  transform  will  apply  to 
the  normal  straight  line  Hough  transform  mentioned  in  Section  H.  References  to 
other  Hough  transform  techniques  will  still  be  preceded  by  the  specific  parameter 
space  identification. 

A.  Accumulator  Technique.  As  mentioned  in  Section  n.D,  the  Hough  transform 
can  be  computed  rapidly  using  an  accumulator  technique.  For  each  edge  point,  a 
sinusoid  was  drawn  in  Hough  space  according  to  Eqn  (2.3) 

p  =  X  cos(0)  +  y  sin(0)  (2.3) 

as  depicted  in  Figure  2.3.  The  Hough  transform  space  with  p  ranging  from  0  to 
181  (128  multiplied  by  the  square  root  of  2)  and  0  ranging  from  0  to  359  was 
created  from  256  x  256  pixel  edge  images. 

The  Hough  transform  was  first  applied  to  simple  shapes  to  best  illustrate 
some  of  its  interesting  qualities.  The  line  in  Figure  4.1  creates  a  local  maxima  at  a 
single  point  in  the  Hough  space  as  shown  in  Figure  4.2.  The  ability  of  tlie  Hough 
transform  to  detect  different  lines  of  the  same  slope  and  also  detect  both  horizontal 
and  vertical  lines  is  displayed  in  Figures  4.3  and  4.4  by  the  Hough  transform  of 
four  line  segments  arranged  in  the  shape  of  a  box.  A  centered  circle  (Figure  4.5) 
creates  maxima  along  the  horizontal  line  with  p  equal  to  the  radius  of  the  circle 
(Figure  4.6)  as  the  Hough  transform  treats  small  arc  segments  of  the  circle  as  line 
segments. 

The  case  involving  the  circle  is  somewhat  interesting  since,  theorectically, 
the  Hough  transform  should  create  a  constant  level  two  out  to  the  radius  of  the 
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Figure  4.3,  Box  Edge  Image 
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Figure  4.4.  Hough  Transform  of  Figure  4.3 

Local  maxima  appear  at  9  =  0®,  90®,  180®  and  270° 
corresponding  to  each  of  the  straight  lines  in  the 
original  image.  Notice  the  wrap  around 
at  9  =  359®  and  9  =  0®. 
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Figure  4.6.  Hough  Transform  of  Figure  4.5 


Hough  transform  of  a  centered  circle  forms  a 
straight  line  with  p  =  circle  radius. 
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The  Hough  transform  was  also  applied  to  the  edge  images  created  from  raw 


An  artifact  of  the  accumulator  technique  is  the  creation  of  many  points  in 
the  Hough  space  with  low  accumulator  values.  This  is  evident  in  the  color 
enhanced  version  of  Figure  4.2  displayed  in  Figure  4.9 


and  its  histogram  presented  in  Figure  4.10. 


Histogram  of  Figure  4.2 
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Figure  4.10.  Histogram  of  Figure  4.2 

A  large  number  of  low  pixel  value  artifacts  are  created  in  the  accumulator  technique. 
For  example,  over  16,000  pixels  have  an  accumulator  value  of  1. 


Hough  transform  images  can  be  thresholded  to  remove  this  artifact  before 
images  are  compared.  It  should  be  noted,  however,  that  useful  information  can  be 
extracted  from  the  artifact  itself.  The  finite  extent  of  this  Hough  space  anifact  can 
help  determine  the  finite  extent  of  the  curve  in  the  edge  image. 


B.  Fourier  Method.  The  Hough  transform,  when  expressed  in  integral  form 
(Eqn  2.8),  is  a  special  case  of  the  two-dimensional  Radon  transform  (Appendix  A) 
were  the  input,  f(x,y),  is  a  binary  function. 


H(0,  p)  =  J J  f(  X,  y)  8(p-  xcos(6)  -  y  sin(0)  ]  dx  dy  (2.8) 


This  integral  relationship  can  be  implemented  through  use  of  Fourier  transforms 
[14:96-100].  The  block  diagram  of  Figure  4.11  outlines  the  procedure  used  in  this 
thesis  to  generate  this  Hough  transform  implementation. 
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Figure  4.11.  Fourier  Transform  Implementation  of  the  Hough  Transform 


A  two-dimensional  Fourier  transform  is  applied  to  the  image  as  defined  in 


Eqn  (4.1). 


F(u,v)  =  JJ  f(x,y) 


-j2'ir(xu  +  yv) 


dx  dy 


(4.2) 


A  Cartesian-to-polar  change  of  coordinates  is  performed  on  the  spatial 
frequency  variables  u  and  v  such  that 


=  yj  and  0  =  tan"^  (  v  /  u ) 


thus, 


_ ,  -  ,  r7  -  j2ir  [  X  r  cos(0)  +  y  r  sin(e)  ] 

F  (  0,  r  )  =  j  j  f(  X,  y  )  e  dx  dy  (4. 


3) 


-OO 


Finally,  by  performing  a  one-dimensional  inverse  Fourier  transform  along 
the  radial  coordinate  of  Eqn  (4.3), 


OO  OO 


H(  0,  p )  «  /  JJf{  X,  y  ) 


-  j  2it  [  X  r  cos(0)  +  y  r  sin(0)] 


dx  dy  e 


j2Trr  p 


dr  (4.4) 


0  -OO 


and  rearranging  terms. 


OO  OO 

Hf  0.  p  )  =  r//  f(  X,  y  )  e 


-j2TTr  [  p  -  xcos(0)  -  y  sin(0) 


dr  dx  dy  (4.5) 


the  inner  integral  simplifies  to  a  Dirac  delta  function  and  Eqn  (2.8)  results. 

Implementing  this  procedure  with  discrete  Fourier  transforms,  the  Hough 
transforms  of  edge  images  in  Figures  4.3,  4.5  and  4.7  are  shown  in  Figures  4.12 
through  4.14. 
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Figure  4.12.  Hough  Transform  of  Figure  4.3  using  Fourier  Method 

Hough  transform  of  the  box  edge  image  in  the  Figtire  4.3  obtained  by 
using  discrete  Fourier  transforms.  Maxima  appear  at  the 
same  positions  as  Figure  4.4. 
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Figure  4.13.  Hough  Transform  of  Figure  4.5  using  Fourier  Method 
Straight  line  in  same  position  as  Figure  4.6. 


Figure  4.14.  Hough  Transform  of  Figure  4.7  using  Fourier  Method 


Notice  that  the  outlines  created  by  the  artifacts  in  Figure  4.7 
are  present  are  also  present  when  Fourier  transforms  are 
used  to  obtain  the  Hough  transform. 

Since  Fourier  transforms  are  invertable,  the  original  edge  images  were 
retrieved  from  the  Hough  transforms  by  using  discrete  Fourier  transforms  and  a 
polar-to-rectangular  coordinate  transformation  to  perform  the  inverse  of  Figure 
4.11. 

The  chapter  outlined  the  distortion  characteristics  of  the  Hough  transform 
due  to  shifts,  rotations  and  translations  of  an  input  object.  Since  these  distortion 
characteristics  follow  well  defined  rules,  they  can  be  used  to  estimate  those  input 
object  parameters  in  the  Hough  space.  There  are  some  interesting  advantages  to 
this  technique  and  they  will  be  discussed  in  the  next  chapter. 


V.  Determining  Scale,  Rotation,  Shift  and  Location 


A..  Distortion  Characteristics.  The  distortion  characteristics  of  the  Hough 
transform  were  described  mathematically  by  Eqns.  (2.4)  through  (2.7). 
Experimental  results  verify  these  equations.  The  effects  of  scales  and  shifts 
(translations)  of  an  input  curve  on  the  Hough  transform  are  best  understood  by 
examining  the  Hough  transform  of  a  circle.  A  centered  circle  of  40  pixel  radius 
and  its  Hough  transform  were  shown  in  Figures  4.5  and  4.6.  The  effect  due  to 
scale  is  evident  when  examining  the  Hough  transform  of  the  circle  of  25  pixel 
radius  (Figure  5.1)  displayed  in  Figure  5.2.  The  effects  of  small  and  large  shifts  of 
the  circle  in  Figure  4.5  (shifts  are  45'  with  respect  to  the  x  axis)  from  center  are 
presented  in  Figures  5.3  through  5.6.  Notice  that  when  the  sinusoidal  shifting  of 
the  p  axis  produces  negative  values,  they  are  shifted  in  theta  by  180  degrees.  This 
accounts  for  band-  of  theta  values  containing  zero  energy  in  Figure  5.6. 


Figure  5.1:  Circle  of  25  pixel  radius 
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Figure  5.2;  Hough  Transform  of  Figure  5.1 


Notice  the  scaling  of  p  when  compared 
to  the  Hough  transform  of  the  40  pixel 
circle  of  Figure  4.6. 


Figure  5.3:  Circle  of  with  Small  Shift 


The  40  pixel  circle  of  Figure  3.3  is  shifted  by  5  pixels. 
The  shift  angle  is  45“  to  the  x  axis. 


Figure  5.6:  Hough  Transform  of  Figure  5.5 

Bands  of  zero  energy,  in  0,  created  by  areas  of  the  Hough 
transform  shifted  by  180“  due  the  creation,  and 
consequent  shifting,  of  negative  p  values. 

B.  Parameter  Estimation  Process.  A  distort  and  compare  routine  was  developed 
to  estimate,  in  the  Hough  transform  domain,  scales,  rotations  and  shifts  of  a 
desired  object  within  an  input  scene.  The  block  diagram  of  Figure  5.7  outlines  the 
process.  Both  input  and  template  Hough  transforms  are  thresholded  to  minimize 
computations.  The  estimator  is  provides  with  search  ranges  and  increment  values 
for  scales,  rotations  and  shifts.  The  template  is  distorted  according  to  Eqn  2.7  for 
each  value  in  the  range.  The  Hough  transform  of  the  input  is  then  compared  to  the 
template.  Scale,  rotation  and  shift  values  corresponding  to  the  best  degree  of 
match  (correlation)  between  input  and  template  are  saved  and  used  to  compute  the 
location  of  the  object  in  the  input  scene. 

Various  methods  can  be  used  to  limit  the  search.  One  is  to  observe  the  zero 
energy  portions  of  the  unthresholded  input  Hough  space  and  omit  shift  values 
within  this  range.  This,  of  course,  is  only  possible  when  the  Hough  space  is 
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determined  by  the  accumulator  technique.  Another  method  is  to  provide  the 
estimator  with  variable  increments.  The  initial  search  could  be  started  with  coarse 
increment  values.  The  increment  values  could  then  be  decreased  as  the  degree  of 
match  (correlation)  increases.  A  final  method  would  be  to  initially  center  the  edge 
images  and  thus  only  search  for  scales  and  rotations  of  the  template. 
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Figure  5.7.  Estimation  Process 

The  thresholded  edge  Hough  transform  image  is 
compared  to  a  template  Hough  transform  image 
that  has  been  distorted  to  simulate  shifts,  rotations 
and  scales  in  the  input  image  space. 
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Table  5.1:  Location  Estimation 

actual 

estimated 

X 

Y 

X 

Y 

5 

15 

5 

15 

24 

25 

24 

25 

36 

-50 

36 

-50 

45 

35 

45 

35 

40 

60 

40 

60 

-60 

60 

-60 

60 

The  Hough  transform  template  for  Table  5.1  developed 
using  a  centered  version  of  the  edge  image  shown  in  Figure  3.3. 

The  edge  image  was  shifted  in  the  scene  by  the  values  in  the 
right  hand  side  of  the  table.  The  Hough  transform  of  the  shifted 
image  estimated  using  the  process  diagramed  in  Figure  5.7. 


The  one  degree  resolution  of  the  Hough  space  was  able  to  determine  location  to 
the  nearest  pixel  in  all  objects  of  reasonable  size.  Location  estimates  were 
developed  by  shifting  a  template  edge  image  in  an  input  scene  and  applying  the 
Hough  transform  to  create  a  new  input  to  the  estimator. 


Table  5.2  Rotation  Estimation 


actual 

estimated 

30 

30 

90 

90 

120 

120 

200 

200 

260 

260 

300 

300 

335 

335 

Table  5.2  shows  Hough  space  estimations  of  the  same  tank  image 
rotated  in  the  input  (edge)  space  by  the  values  shown. 


Table  5.3c  Scale  Estimation 


tank  1 

actual  1 

estimated 

d3033  1 

1.0 

1.0 

d3108 

0.7 

0.7 

d3042 

0.4 

0.4 

d3074 

0.3 

0.2 

Rotation  estimates  were  also  within  the  increment  resolution  of  the  search.  The 
accuracy  of  the  location  and  rotation  estimates  indicates  that  a  Hough  space  of  less 
resolution  would  be  sufficient  to  characterize  most  input  scenes.  The  scale 
estimates  used  similar  tank  edges  of  differing  sizes.  Actual  size  of  these  cases 
were  determined  by  length  and  width  comparisons  with  the  template.  The  scale 
estimator  was  unable  to  accurately  determine  scales  below  approximately  0.4  due 
to  collapse  of  the  Hough  space  about  p  =  0. 

This  chapter  has  demonstrated  how  the  distortion  characteristics  of  the 
Hough  transform  due  to  rotations,  scales  and  shifts  of  an  input  object  can  be  used 
to  estimate  these  parameters  in  the  Hough  space.  The  following  chapter  contains 
an  overall  discussion  of  methods  used  in  this  thesis  and  there  relative  performance. 


VI.  Discussion 


A.  Accumulator  Method.  The  accumulator  method  provided  an  efficient  means  of 
producing  the  Hough  transform.  The  low-level  noise  artifacts  present  a  major 
limitation  to  employing  the  transform  to  more  than  one  segmented  edge  object  at  a 
time.  However,  by  applying  the  transform  to  each  area  and  thresholding  them 
separately  to  remove  the  artifacts,  the  method  can  be  used  to  characterize  differing 
image  outlines.  In  well  segmented  images,  the  artifacts  actually  proved  to  be 
useful  in  reducing  the  range  of  the  search  in  the  parameter  estimation  process. 
Segmentation  noise  has  a  non-linear  effect  on  the  Hough  transform.  Stray  pi.xels 
present  in  the  outer  portions  of  the  image  produce  more  noise  in  the  transform 
than  those  toward  the  center  since  the  sinusoid  it  produces  is  weighted  by  the 
distance  of  the  pixel  from  the  origin.  Thus,  adequate  median  filtering  is  essential 
in  the  initial  segmentation  process  before  application  of  the  Hough  transform. 

B.  Parameter  Estimation.  When  artifact  information  can  be  employed,  or  when 
the  input  image  is  centered  initially  to  narrow  the  range  of  the  search,  the 
parameter  estimation  method  can  distinguish  between  differing  images  with 
reasonable  efficiency.  Size  and  rotation  information  are  also  extractable  in  the 
estimation  process.  Rotation  increments  near  6  degrees  and  size  increments  of 
approximately  0.1  proved  quite  adequate  for  estimating  locations  of  a  template 
among  multiple  objects  in  an  image.  Rotation  and  shift  were  determined  e.xactly 
(or  to  the  resolution  of  the  search  increment)  while  scale  estimates  proved  accurate 
when  the  input  image  was  at  least  0.4  of  the  template. 


33 


A, 


V  ••  Si 

•  V  V 

'•*.  ••  s 

•  V/v.V 


'■nA’/.’A 


V  V  I 

V.V  ■«*"*" 


C.  Fourier  Method.  Use  of  discrete  Fourier  transforms  and  a  coordinate 
transformation  to  produce  the  Hough  transform  correlated  highly  with  the 
accumulator  method.  The  ’butterfly’  patterns  of  line  segments,  as  well  as  the  more 
intricate  variations  in  complex  shapes  matched  those  produced  by  the  accumulator 
method.  The  Hough  transform  of  a  circular  object  was  highlighted  in  Section  IV.A 
noting  the  maximum  created  in  the  Hough  space  at  p  equal  to  the  radius  of  the 
circle  when  the  accumulator  method  was  used.  Interestingly,  this  maxima  is  also 
present  in  the  Fourier  created  Hough  space.  Thus,  the  Fourier  method  treats  small 
arc  segments  in  the  same  fashion  as  the  accumulator  method.  Notice  that  if  a 
two-dimensional  delta  function  is  inserted  into  Eqn  (2.5),  or  Eqn  (4.8),  the  result 
is  a  line  mass  distributed  along  a  sinusoid.  In  both  the  continuous  and  discrete 
(accumulator)  implementations,  each  point  in  the  input  contributes  to  a  sinusoid  in 
the  Hough  space.  Therefore,  what  is  commonly  referred  t&  as  the  accumulator 
method  should,  in  fact,  be  called  the  discrete  Hough  transform.  The  next  chapter 
outlines  some  attempts  to  achieve  the  continuous  representation  optically  which, 
unfortunately,  proved  unsuccessful. 


VII.  Suggested  Optical  Implementation. 


A.  Flawed  Approach.  The  effort  to  improve  upon  current  optical 
implementation  methods  produced  no  useful  results.  The  major  obstacle  to  a 
continuous  real-time  optical  reproduction  of  the  Hough  transform  is  the 
rectangular-to-polar  coordinate  transformation  required  on  the  complex  field  in 
the  image  Fourier  transform  plane.  One  method  considered  made  use  of 
computer-generated  interferograms  (Appendix  C)  to  implement  the  coordinate 
transformation.  The  computer-generated  interferograms,  unfortunately,  are 
effective  only  on  planar  wavefronts. 

B.  New  Information.  A  recent  article  by  Jensen  et  al  [18]  suggests  a  variation 
on  the  CGH  method.  Jensen  recorded  the  phase  function  of  the  [x,  y]  -  to  -  [6. 
lnV(x2  +  y2)]  Fourier  transform  coordinate  transform  (see  Appendix  B)  on 
dichromated  gelatin.  With  their  phase  filter,  they  effected  the  coordinate 
transformation  on  a  complex  wavefront.  Jensen  et  al  claims  that  derivation  of  Eqn 
(B.8)  will  hold  for  the  Fourier  transform  of  a  real  input  and  has  published  data  on 
performing  the  coordinate  transformation  in  the  Fourier  transform  plane.  They 
state  that  as  long  as  the  variations  in  the  input  are  small  compared  to  the  variations 
in  the  phase  filter,  the  coordinate  transformation  can  be  performed.  If  this  is 
indeed  the  case,  the  following  optical  scheme  could  be  implemented, 

C.  Suggested  Approach.  By  substituting  Jenson  phase  filters,  or  possibly  one-bit 
phase  filters  using  chemically  etched  glass,  for  the  computer-generated  filters,  a 
possible  approach  for  generating  a  continuous  Hough  (Radon)  transform  is 
suggested.  In  the  implementation,  shown  in  Figure  7.1,  the  Fourier  transform  of 
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the  input  is  multiplied  by  the  phase  function  of  Eqn  (B.8)  in  plane  P2-  This 
wavefront  is  Fourier  transformed  to  produce  the  [0,  In  r]  transformation.  This 
wavefront  is  multiplied  by  the  one-dimensional  exponential  phase  filter 


4),(0,  0=  exp(l)  (7.1) 

This  is  again  Fourier  transformed  to  produce  the  exponential  shifting  of  the  In  r 
coordinate.  Finally,  a  cylindrical  lens  is  suggested  to  do  the  final  one-dimensional 
Fourier  transformation. 


input 
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phase 
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conjugate 
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filter,  Eqn  (7.3),  is  place  in  P^  and 


the  phase  filter 


4)1  (0,  r)  =  ln(r) 


(7.3) 


4)(0,  1)  =  exp(l)  sin(0)  (7.4) 

used  to  produce  the  [0,  l]-to-[x,  y]  coordinate  transformation.  This  phase  function 
satisfies  the  four  partial  differential  equations  of  Appendix  C. 


lens  (Appendix  B). 

This  optical  scheme  is  presented  as  a  suggestion  for  further  research.  Its 
feasibility  depends  highly  on  the  accuracy  of  the  Jensen  et  al  development.  This 
point  is  reiterated  in  the  conclusion  to  this  thesis  which  follows  in  next  chapter. 
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VIII.  Conclusion 


This  thesis  has  tested  the  feature  extraction  properties  of  the  Hough  trans¬ 
form  on  simple  images  and  real  data.  Experimental  results  have  demonstrated 
that  Hough  transform  techniques  to  determine  shifts,  rotations  and  scales  of  an 
input  object  in  the  Hough  space  can  be  easily  implemented  through  a  distort-and- 
compare  process.  The  accuracy  in  determining  these  distortion  parameters  were 
primarily  influenced  by  the  resolution  of  the  search  increments  and  of  the  gener¬ 
ated  Hough  space.  When  the  accumulator  technique  is  used,  information  to  re¬ 
duce  the  range  of  the  search  can  be  extracted  from  low  pixel  value  artifacts  in  the 
Hough  transform  of  an  input  edge  image. 

A  continuous  implementation  of  the  Hough  transform  using  Fourier  trans¬ 
forms  was  implemented  successfully.  This  points  favorably  to  a  possible  continu¬ 
ous  optical  implementation.  The  main  obstacle  is  producing  the  necessary  rectan- 
gular-to-polar  coordinate  transformation. 

Attempts  to  produce  a  continuous  optical  implementation  of  the  Hough  trans¬ 
form  using  computer-generated  interferogram  coordinate  transformations  proved 
unsuccessful.  In  order  for  this  approach  to  succeed,  the  complex  field  requiring 
coordinate  transformation  must  be  separated  and  processed  individually.  A  recent 
article  by  Jensen  et  al  may  provide  an  answer  to  this  problem  of  a  continuous 
implementation.  The  assumptions  made  in  his  coordinate  transformation  deriva¬ 
tion  (Apfjendix  B)  may  be  too  restrictive  for  this  application.  Thus,  it  is  only 
presented  as  a  suggestion  for  further  research. 


Appendix  A:  Radon  Transform 


This  appendix  is  limited  to  discussion  of  the  two-dimensional  Radon 
transform.  An  extension  to  higher  orders  can  be  found  in  Deans  [14]. 

1.  General  Definition.  The  Radon  transform  is  formed  by  mapping  the 
projection  of  some  function  along  all  possible  lines  in  a  plane.  Thus,  the  line 
integral  of  Eqn  (A.l)  describes  the  general  Radon  transform. 

R(e,  p)  =  f  f(  X,  y)  dl  (A.l) 

•'l 

Thus,  the  first  restriction  is  that  this  integral  must  exist.  In  his  paper, 
published  in  1917  [14:204-21 7]  Radon  showed  that  the  transform  must  exist  if  the 
following  conditions  are  meet: 


a)  f(x,y)  is  continuous 

b)  the  integral  of  Eqn  (A.2)  converges 


/ 

-  00 


f(  x,  y  ) 

Vrx2+  y2) 


dx  dy 


c)  and,  for  in  the  plane  the  limit  in  Eqn  (A.3)  holds 


(A.2) 


Itt 

lim  r  f(  X  +  r  cos(0),  y  +  r  sin(0)  )  d0  (.A.3) 

r  oo  2-7r  Q 


These  conditions  can  be  satisfied  for  most  functions.  For  instance,  if  the 
function  is  bounded,  the  last  two  restrictions  will  be  satisfied.  With  discontinuous 


function,  the  remaining  condition  can  still  be  satisfied  by  integrating  over 
continuous  segments.  Starting  with  the  parameterization  of  Eqn  (2.3),  and 
superimposing  a  new  coordinate  system  (p,  s)  over  Figure  2.2  where  L  is 
perpendicular  to  the  p  axis  (Figure  A.l), 


Figure  A.l.  Line  Through  f(x,  y)  [14:57] 


each  point  along  L  is  determined  by 


X  a  p  cos(0)  -  s  sin(0) 
y  =  p  sin(0)  +  s  cos(0) 


Thus.  Eqn  (A.l)  is  explicitly  defined  as 


M  -  * 


R(  0,  p  )  *  f  f(  P  cos(0)  -  s  sin(0),  p  sin(0)  +  s  cos(0)  )  ds  (A.5) 

-  OO 


Alternatively,  the  Dirac  delta  function  can  be  employed  to  select  each  possible  line 
L  of  Figure  A.l  and  the  more  familiar  Eqn  (A.6)  results. 


OO 

R(0,  p)  =  rf(x,y)8[p-x  cos(0)  -  y  sin(0)  ]  dx  dy 


(A.6) 


Appendix  B:  Coordinate  Transformations 


Information  in  this  appendix  was  extracted  from  Casasent-Psaltis  [16], 
Casasent  et  al  [15],  Bom-Wolf  [22],  Jenson  et  al  [19]  and  Bryngdahl  [18]. 

1.  Fourier  Transform  Approach.  Casasent  and  Psaltis  [16]  introduced  the 
concept  of  Fourier  transform  coordinate  transformations.  The  approach  states  that 
an  input  function,  f(x,  y),  can  be  multiplied  by  an  appropriately  selected  phase 
filter  function  to  produce  an  output  function  whose  intensity  is  ’proportional’  to 
some  desired  f(u,  v)  where  u  and  v  are  functions  of  x  and  y.  This  is  stated 
mathematically  in  Eqn  (B.l)  and  optically  implemented  in  Figure  B.  1. 


f  f  f(  X,  y  )  exp{  j  ^  y  exp{  -j 


211  (xu  +  yv) 


}  dx  dy 


(B.l) 
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Figure  B.l.  Optical  Implementation  of  Coordinate  Transform  Filter 


In  order  to  determine  the  limitations  on  choosing  a  4)(  x,  y  ),  a  solution  to 
Eqn  (B.l)  must  be  derived.  Methods  for  solving  this  double  integral  are  provided 
by  Born-Wolf  [22]  and  Jenson  et  al  [16].  Although  outlined  quite  differently,  both 
methods  follow  the  same  reasoning.  Therefore,  only  the  method  published  in 
Born-Wolf  [22]  will  be  broached. 

The  method  discussed  in  Appendix  HI  of  Bom  and  Wolf  [22]  used  an 
as>  mtotic  approximation  to  solve  double  integrals  of  the  form 


OO 

f  f(  X,  y  )  exp  {  j  kh(  X,  y  )  }  dx  dy  (B.2) 

-OO 

This  approach  limits  h(  x,  y  )  to  a  real  function.  In  general  terms,  the  approach 
holds  that  ’’contributions  to  the  asymtotic  expansion  (of  h(  x,  v  ))  come  only  from 
regions  in  the  vicinity  of  certain  critical  points.”  Three  types  oi  v  ritical  points’  are 
discussed;  however,  only  the  case  called  ’a  critical  point  of  the  iirst  kind’  pertains 
to  this  application.  At  each  critical  point  where, 


SJl  =  SJl  =  0 

8  X  8  y 

provide  a  necessary  condition  for  choosing  4)(  x,  y  ).  Since 


(B.3) 


One-dimensional  image 


3.  One-Dimensional  Image  Modification. 
transformation  filters  have  no  requirement  such  as  Eqn  (B  5).  Optical 
implementation  of  the  one-dimensional  filters  is  identical  to  two-dimensional 
coordinate  transformations.  The  exponential  image  modification  filter  is  shown  in 
Figure  B.3  using  2ir  contour  lines. 


Figure  B,3.  Exponential  Image  Modification  Phase  Filter 

One-dimensional  phase  filter  for  exponential 
expansion  of  input  wavefront.  Plot 
generated  from  2tt  contours  of 
the  phase  htnction. 


Appendix  C:  Computer-Generated  Interferograms 

Information  in  this  appendix  was  extracted  from  Lee  [17:152-154]  and 
Casasent  et  al  [15]. 

1 .  Definition.  The  computer-generated  interferogram  is  an  optical  technique  for 
producing  an  arbitrary  complex  field  using  a  binary  amplitude  transmittance.  The 
complex  field  is  created  on  a  computer  and  an  amplitude  transmittance  produced 
by  computing  the  interference  pattern  of  a  reference  wave  and  the  complex  field. 
The  result  is  then  output  to  a  plotting  device  and  photo-reduced  onto  an  optical 
slide  for  use. 

2.  Phase-Only  Application.  Computer-generated  interferograms  are  normally 
used  to  produce  phase-only  filters.  In  this  application,  the  interference  pattern  of 
a  desired  phase  function  and  an  off-axis  uniform  amplitude  plain  reference  wave 
is  computed  according  to  Eqn  (C.l). 

t(x,y)  =  0.5  { 1  +  cos[  217  a  X  -  4>(  X,  y  )  ]  }  (C.l) 

Maxima  in  this  equation,  which  occur  when 

2tt  ot  \  -  4)(x,y)  =  217  n  ,  n  =  0,1,2,... 

form  the  binary  amplitude  transmittance.  Output  plots  can  then  be  photo-reduced 
onto  the  appropriate  optical  media.  When  this  transmittance  is  illuminated  with  a 
planar  wavefront,  the  desired  combination  of  phase  filter  and  input  wavefront 
appear  off-axis  at  an  angle  proportional  to  that  of  the  reference  wave  Ck:.  Figure 


Appendix  D:  Thesis  Resources 


This  appendix  lists  the  primary  resources  used  to  complete  this  thesis  effort. 
The  doppler  data,  computer  resources  and  graphic  resources  were  all  made  avail¬ 
able  by  previous  AFTT  research  sponsors. 


•A.  Doppler  Data.  Doppler  images  were  supplied  by  the  Air  Force  Wright  Aero¬ 
nautical  Laboratory,  Avionics  Laboratory  (AFWAL/ALTC). 


B.  Computational  Resources.  The  bulk  of  the  computer  processing  was  per¬ 
formed  on  the  AFTT  Information  Sciences  Laboratory  (ISL)  VAX  11/780  computer. 
Contouring  for  generation  of  computer-generated  interferograms  was  performed 
on  an  AFTT  ICC  Elxsi  computer  and  AFTT  Physics  Department  Sun  n  Workstation. 


C.  Graphic  Resources.  Edge  and  Hough  transform  images  were  displayed  on  the 
.AFIT  Department  of  Electrical  and  Computer  Engineering  Evans  and  Sutherland 
PS340  Graphics  Workstation.  Computer-generated  interferograms  were  output  on 
the  .-^FIT  Physics  Department  Imagen  laser  printer  and  photo-reduced  in  the 
photo-processing  laboratory  located  in  the  AFTT  Headquarters  building. 


Appendix  E:  VAX  Programs 


The  Ada  code  developed  during  this  thesis  effort  which  specifically  apply  to 
the  material  presented  in  this  document  are  listed  in  this  appendix.  Several  pack¬ 
ages  used  in  some  of  the  specifications  but  not  list  are  edited  versions  of  code 
developed  in  previous  thesis  work.  They  include  the  following; 


Package  Name 

Original  Name 

Original  Author 

Data_Standard 

Data_Standard 

Tong  [20] 

Data_Conversions 

Data_Conversions 

Tong  [20] 

FileJO 

FileJO 

Tong  [20] 

SSIJO 

FileJO 

Ruck  [23] 

Fourier2 

Fourier_Transform_Handler 

Kobel-Marrin  [22 

The  .Ada  programs  developed  for  this  thesis  are  listed  in  the  following  pages  of  this 
appendix. 


PACKAGE  HFl 


—  Purpose:  Digital  filters  for  straight  line  Hough  images 


with  DATA_STANDARD:  use  DATA_STANDARD ; 

with  floatjnath_lib;  use  float_math_lib; 

with  text_io;  use  text_io; 

with  integer_text_io ;  use  integer_text_io ; 

package  HFl  is 

procedure  COMPARE (  Inputl:  in  DATA_STANDARD.Image_ArraY_2d; 

Input2;  in  DATA_STANDARD. Image_Array_2d; 
Match  ;  out  integer  ) ; 


procedure  DISTORT (  Input  :  in  DATA_STANDARD . Image_ArraY_2d ; 

Rotate ;  in  integer : 

Radius;  in  integer; 

Angle  :  in  integer; 

Scale  :  in  float; 

Output:  in  out  DATA_STANDARD.Image_ArraY_2d  ); 


procedure  ESTIMATE (  Input 

Template 

Rot_.Min 

Rot_Max 

Rot_Delta 

IRad_MIn 

Rad_Max 

Rad_Delta 

Ang_Min 

Ang_Max 

Ang_Delta 

Scale_Min 

Scale_Max 

SccLle_Delta 

X_Location 

y_Locatian 

Best_Rotate 

Best_Scale 

Best  Match 


in  DATA_STANDARD . Image_ArraY_2d ; 

in  DATA_STANDARD .  Iinage_ArraY_2d  ; 

in  integer; 

in  integer; 

in  integer; 

in  integer; 

in  integer : 

in  integer; 

in  integer; 

in  integer; 

in  integer : 

in  float; 

in  float; 

in  float ; 

out  integer; 

out  integer; 

out  integer; 

out  float; 

out  integer  ) ; 


—  FUNCTION:  Coi^are  tvro  iamges 

—  INPUTS  ( 1 )  Ii^tl  -  Fi2?st  Image 

( 2 )  Ir^t2  -  Second  Image 

—  OUTPUT  ( 1 )  Match  -  Percent  match  { integer  0  -  95 ) 

_ m******** ****************  ***************************************************** 

procediire  COMPARE (  Input 1:  in  DATA_STANDARD . Image_Array_2d ; 

Input  2 :  in  DATA_STANDARD . Image_Array_2d : 

Match  :  out  integer  )  is 


Add,  Sum 
Count 

SIZE_ERROR 
begin  —  COMPARE 


float  :=  0.0; 
integer  :=  1; 
exception; 


—  [1]  Check  for  size  differences 

if  Inputl 'first(l)  /*  Input2’first(l)  or  Ir^tl ' last ( 1 )  /-  Input2 ' last ( 1 )  then 
raise  SIZE_ERROR; 

elsif  Inputl 'first{2)  /=  Ii^t2 '  f irst ( 2 )  or  Iiputl'last(2)/=  Input 2 '  last ( 2 )  then 
raise  SIZE_EF?ROR; 
end  if; 

—  [2]  Conpare  images 

for  Theta  in  Inputl 'first(l) .  .Ir^tl 'last{l) 
loop 

for  Rho  in  Ir5Jutl'first(2) .  .Ir^tl'last(2) 
loop 

if  Inputl (Theta, Rho)  >  0  then 

if  Input2( Theta, Rho)  >  0  then 

if  Inputl (Theta, Rho)  <  Ir¥Mt2( Theta, Rho)  then 

Add  :=  float(  Ir^jutl{ Theta, Rho)  )  /  float(  Input2( Theta, Rho)  ); 
else 

Add  :=  float(  Input2( Theta, Rho)  )  /  float(  Input 2 (Theta, Rho)  ); 
end  if; 

Sum  :=  Sum  +  Add; 
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end  if; 


Count  ; =  Count  +  1 ; 


end  if; 
end  loop; 
end  loop; 


if  Count  >  1  then 

Count  : =  Count  -  1 ; 
end  if; 


Match  integer (  100.0  ♦  Sum  /  float (Count)  ); 


exception 

vghen  SIZE_ERROR  => 
new_line; 

put_line("  *•*  SIZE  ERROR  **•"); 
put_line("  Images  of  unequal  sizes."); 


erxi  COMPARE; 

_ ^:^m**m********* *************************************************************** 


DISTORT 


FUNCTION:  Adjust  Hough  space  for  a  give  shift  in  tlie  image  space 
INPUTS: 


(1)  Ir^t  -  Input  straight  line  Hough  image 

(2)  Rotate  -  Rotation  in  image  space  (degrees) 

(3)  Radius  -  Shift  radius  in  image  space  (pixels) 

(4)  Angle  -  Shift  angle  in  im^e  space  (degrees) 

(5)  Scale  -  Scale  in  image  space  (fraction) 


—  OUTPUTS:  (1)  Output  -  Adjusted  straight  line  Hough  image 

_ ^mm*************************************************** ************************ 


procedure  DISTORT ( 

Input  : 

in 

Rotate : 

in 

Radius: 

in 

Angle  : 

in 

Scale  : 

in 

Output: 

in  out 

Rad,  Ang 

:  float ; 

Rho_In,  Theta_In 

;  float ; 

Rho_0ut ,  Theta_0ut 

:  integer; 

begin  —  DISTORT 


54 


—  [ 1 ]  Set  Output  to  zero 

for  Theta  in  Output ' first ( 1 ) . . Output ' last ( 1 ) 
loop 

for  Rho  in  Output ' first (2) . . Output ' last ( 2 ) 
loop 

Output (Theta, Rho)  ;=  0; 
end  loop; 
end  loop; 

—  [2]  Distort  image 

Rad  :=  float (  Radius  ); 

Ang  :=  float (  Angle  ); 

for  Theta  in  Input ' f irst( 1) . .Input 'last (1) 
loop 

foi-  Rho  in  Ir^t '  first  ( 2 )..  Input '  last  ( 2 ) 
loop 

if  Input (Theta, Rho)  >  0  then 

Rho_In  :=  float (Rho); 

Theta_In  :=  float ( Theta ) ; 

Rho_Out  :=  integer(  (Rho_In  +  Rad  *  cosd ( Theta_In  -  Ang))  •  Scale  ) 

if  Rho_0ut  <  0  then 

Iheta_0ut  :=  Theta  -  Rotate  +  180; 

Rho_0ut  :=  -Rho_0ut; 
else 

Theta_0ut  Theta  -  Rotate; 
end  if; 

if  Theta_Out  >  359  then 

Theta_0ut  :=  Theta_0ut  -  360; 
elsif  Theta_0ut  <  0  then 

Theta_0ut  :=  'lheta_0ut  +  360; 
end  if; 

if  Input (Theta, Rho)  >  Output (Theta, Rho)  then 
Output ( Theta_0ut , Rho_0ut )  ; =  Input ( Theta , Rho ) ; 
end  if; 

end  if; 
end  loop; 
end  loop; 


end  DISTORT; 


ESTIMATE 


—  FUNCTION:  Estimate  location  and  scale  of  tenplate  in  input  image 

_ innL*^*mm**********m****»*******’*********************************************** 


procedure  ESTIMATE ( 


Input 

in 

DATA_STANDARD 

Template 

in 

DATA_STANDARD 

Rot_Min 

in 

integer ; 

Rot_^fex 

in 

integer; 

Rot_Delta 

in 

integer; 

Rad_MIn 

in 

integer; 

Rad_Max 

in 

integer; 

Rad_Delta 

in 

integer ; 

Ang_Min 

in 

integer; 

Ang_Max 

in 

integer; 

Ang_Delta 

in 

integer; 

Scale_Min 

in 

float; 

Scale_Max 

in 

float; 

Scale_Delta 

in 

float; 

X_Location 

out 

integer ; 

Y_Location 

out 

integer; 

Best__Rotate 

out 

integer; 

Best_Scale 

out 

float; 

Best__Match 

out 

integer  )  is 

Output 

Done 

Rotate 

P^adixis 

Angle 

Scale 

Match 

Best_Rad 

Best_Ang 

Best  Mat 


DATA_STANDARD .  Image_Array_2d  (  Ir^t '  first  ( 1 ) . .  Input '  last  ( 1 ) , 

Input' first! 2) . .Input 'last(2)  ) 
DATA_STANDARD.Bit_Array(  1..4  )  :=  (  true,  true,  true,  true  ) 


integer 

integer 

integer 

float 

integer 

integer 

integer 

integer 


=  Rot^Min; 

=  Rad__Min: 

=  Ang__Min; 

=  Scale_Min; 
=  0; 

=  0 
=  0 
=  0 


begin  —  ESTIMATE 

BestScale  :=  Scale; 

while  Rotate  <=  Rot_Max 
loop 

while  Radius  <=  Rad_Max 
loop 


while  Angle  •<=  Ang_Max 


loop 


while  Scale  <=  Scale_Max 
loop 


DISTORT (  Tenplate,  Rotate,  Radi  vis.  Angle,  Scale,  Output  ) ; 
COMPARE (  Output,  Input,  Match  ); 

if  Match  >  Best  Mat  then 


Best_Mat 
Best_Scale 
Best_Rad 
Best_Ang 
Best_Rotate 
end  if; 


=  Match; 

=  Scale; 

=  Radius; 
=  Angle; 

=  Rotate; 


Scale  ;=  Scale  +  Scale_Delta; 
end  loop; 

Scale  :=  Scale_Min; 

Angle  :=  Angle  +  Ang_Delta; 


"t 
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end  loop; 

Angle  :=  Ang_Min; 

Radius  :=  Radius  +  Rad_Delta; 

end  loop; 

Radius  :=  Rad_Min; 

Rotate  :=  Rotate  +  Rot_Delta; 


end  loop; 


Best_Match 
X_Location 
Y  Location 


=  Best_Mat; 

=  integer (  float (Best_Rad)  ♦  cosd(  float (Best_Ang)  )  ); 
=  integer(  float (Best_Rad)  *  sind(  float (Best_Ang)  )  ) ; 


end  ESTIMATE; 
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— ****************************  gjjjj  package  HFl  ********************************* 
end  HFl ; 


•  I 


with  DATA_STANDARD; 
with  FILE_IO; 
with  SSI_IO; 
with  HFl; 

with  integer_text_io; 
with  float_text_io; 
with  text_io; 

procedure  ESTIMATOR  is 


use  DATA_STANDARD; 
\ise  FILE_IO: 
use  SSI_IO; 
xise  HFl; 

use  integer_text_io ; 
\ise  f loat_text_io ; 
use  text  io; 


Input ,  Template 
Input_Count ,  TanpIate_Count 
Input_Name ,  Template_Naiue 

Choice 

Rot_Min 

Rot_Max 

Rot_Delta 

Rad_Min 

Rad_Max 

Rad_Del -a 

Ang_Min 

Ang_Max 

Ang_Delta 

Scale_Min 

Scale_Max 

Scale_Delta 

X_Location 

Y_Location 

Best_Rotate 

Best_Scale 

Best_Match 

begin  —  ESTIMATOR 

—  [1]  Get  input  HS  file 


new_line; 

put_line("  ENTER  NAME  OF  INPUT  HS  SSI  FILE"): 
GET_FILENAME(  Input_Name,  Input_Count ) ; 
READ_SSI_IMAGE(  Input_Name,  Input_Count,  Input) 


DATA_STANnARD.Image_Array_2d(0. .359,  0. ,181) 
natural ; 

DATA_STANDARD . Line_Fonn ; 
character ; 


integer 

integer 

integer 

integer 

integer 

integer 

integer 

integer 

integer 

float 

float 

float 

integer 

integer 

integer 

float 

integer 


—  [2]  Get  Template 


new_line; 

put_line(''  ENTER  NAME  OF  TEMPLATE"); 

GET_FILENAME(  TaTiplate_Name ,  Template_Count  ); 

READ_SSI_IMAGE(  Template_Name ,  Tenplate_Count ,  Tenplate  ) ; 

—  [3]  Get  Range 

GET_LOOP: 

loop 

new_line(3) ; 

put_line("  *•*  HOUGH  SPACE  ESTIMATOR  ***  "); 
new_line; 

put("  Ir^t  Name  =>  ");  put ( Input_Name ( 1 .  . Ir^t_Ccfunt ) ) ; 

new_line; 

put{"  Taiplate  ^feune  =>  ");  put{TempIate_Name(l .  .Tanplate_Count) ) 

new_line(2) ; 

put("  [a]  Min  Rotation  Value.  =  ");  put(Rot_Min,3, 10) ;  new_line; 

put("  [bj  Max  Rotation  Value  =  ");  put(Rot_Max,3, 10) ;  new_line; 

put("  [c]  Rotation  Delta  =  ");  put{Rot_Delta,3, 10) ;  new_line(2); 

put("  [d]  Min  Radius  Value  =  ");  put(Rad_Min,3, 10) ;  nevg'_line; 

put("  [e]  Max  Raditis  Value  =  ");  put{Rad_Max,3,10) ;  new_line; 

put("  [f]  Radius  Delta  =  ");  put(Rad_Delta,3, 10) ;  new_line(2) ; 

put("  [g]  Min  Angle  Valvie  =  "):  put(Ang_Min,3, 10) ;  new_line; 

put("  [h]  Max  Angle  Value  =  ");  put(Ang_Max,3, 10) ;  new_line; 

put("  [i]  Angle  Delta  =  ");  put(Ang_Delta,3,10) ;  new_line(2); 

put("  [j]  Min  Scale  =  ");  put(Scale_Min, 1 , 2, 0) ;  new_line; 

put("  [k]  Max  Scale  =  ") ;  put(Scale_Max,l,2,0) ;  new_line; 

put("  [1]  Scale  Delta  =  ");  put(Scale_Delta, 1 , 2,0) ;  new_line(2); 

put_line( "  [z]  Perform  Estimation. . . ") ; 
put("  ENTER  CHOICE  >  ");  get ( Choice ) ;  new_line; 

case  Choice  is 

when  'a'  (  'A'  => 

put("  Enter  min  rotation  value  >  "); 
get(Rot_Min) ; 

vrfien  'b'  I  'B'  => 

put("  Enter  max  rotation  value  >  "); 
get{Rot_Max)  ; 

when  'c'  |  'C'  => 

put("  Enter  rotation  delta  >  "); 
get (Rot  Delta); 


when  'd'  (  'D'  => 

put ( "  Enter  min  radius  value  >  " ) ; 
get(Rad_Min)  : 

when  ' e '  (  ' E '  => 

put ( "  Enter  max  radius  value  >  " ) ; 
get(Rad_Max) : 

when  ' f  I  ' F '  => 

put  ( "  Enter  radixis  delta  >  " ) ; 
get ( Rad_Del ta ) : 

when  'g'  |  'G'  => 

put ( "  Enter  min  angle  value  >  " ) ; 
get(Ang_Min) ; 

when  'h'  [  'H'  => 

put{"  Enter  max  angle  value  >  "); 
get(Ang_Max)  ; 

when  ' i '  |  ' I '  => 

put("  Enter  angle  delta  >  "); 
get(Ang_Delta) ; 

when  ' j '  I  ' J '  => 

put  ( "  Enter  min  scale  value  >  " )  ; 
get ( Scale_Min) ; 

when  'k'  |  'K'  => 

put ( "  Enter  max  scale  value  >  " ) ; 
Sfet{Scale_Max) ; 

vAien  '  1 '  I  '  L '  => 

put{”  Enter  scale  delta  >  "); 
get { Scale_Del ta ) ; 

when  'z'  |  'Z'  => 
exit  GET_LOOP; 

when  others  =>  null; 

end  case; 

end  loop  GET_LOOP; 

—  [4]  Perform  estimation 

put_line( "  asREM  -  Estimating. . .  " )  ; 

ESTIMATE {  Input,  Template, 

Rot_Min,  Rot_Max,  Rot_Delta 

Rad  Min,  Rad  Max,  Rad  Delta 


Ang_Min,  Ang_Max,  Ang_Delta, 

Scale_Min,  Scale_Max,  Scale_Delta, 

X_Location,  Y_Location,  Best_Rotate,  Best_ScaIe,  Best_Match  ) 

—  [5]  Print  output 

new_Iine(2) ; 
put("  X 
put("  Y 
put ( "  Rotation 
put { "  Scale 
put ( "  Match 

new_line(2) ; 

put_line("  *•*  END  ESTIMATOR  ***"); 
new  line(2) ; 


=  ");  put(X_Location,  3,  10);  new_line: 

=  ");  put(Y_Location,  3,  10);  new_line; 

=  "):  put(Best_Rotate.  3,  10);  new_Iine; 

=  ");  put(Best_Scale,  1,  2,  0);  new_line; 
=  ");  put(Best_Match,  3,  10);  new_line; 


end  ESTIMATOR; 


FOURIER  RADON 


—  Purpose:  To  generate  a  Radon  Transform  through  a  2D  Fourier  Transform 

—  followed  by  a  inverse  Fourier  Transform  along  the  radial  direction. 

_ if:>Ht:tm*****************************************m********************  *********** 

with  DATA_STANDARD;  use  DATA_STANDARD ; 

with  DATA_CONVERSIONS;  use  DATA_CONVERSIONS : 

with  F0URIER2;  use  F00RIER2; 

with  SSI_IO;  use  SSI_IO; 

with  FILE_IO;  use  FILE_IO; 

with  text_io;  use  text_io; 

with  integer_text_io ;  use  integer_text_io ; 

with  float_math_lib;  use  float_math_lib; 

proced\are  FOURIER_RADON  is 

Input,  FT_Output 
Output,  FT_RT_Output 
In_Work 
Out_Work 
Rho_Wbrk 
Output_Histogram 
Out^t_Histogram_float 
Name 

Direction 

X_Center,  Y_Center, 

Nearest_X,  Nearest_Y 
Max_Output 
Min_Output 
Temp 

Tejnp_Out , 

X_Arg,  Y_Arg, 

X_Arg_Cos,  X_Arg_Sin, 

Y_Arg_Cos ,  Y_Arg_S in , 

In_fteal,  In_ljng, 

Real_Shift,  Img_Shift 
Char_Count 
Answer 
Edge_Value 

begin  ~  FOURIER_RADON 

—  [1]  Get  input  image 
new_line; 

put_line("  ENTER  SSI  INPUT  IMAGE"); 

GET_FILENAME(  Name,  Chau:_Count  )  ; 


DATA_STANDARD . Image_Array_2d (  1 . . 256 ,  1 . . 256  ) ; 
DATA_STANDARD . Image_Array_2d (  0 . . 359 ,  0 . . 18 1  ) ; 
nATA_STANDARD . Ccn?)lex_Array_R2d {  1..256,  1..256  ); 
DATA_STANDARD . Cai¥)lex_AiTay_R2d (  0..359,  0..181  ); 
DATA_STANDARD .  Coinplex_Array_Rld  {  0 .  ,  255  ) ; 
DATA_STANDARD.Image_Array_ld(  0..255  ); 
DATA_STANDAFD.Array_ld(  0..255  ); 

DATA_STANDARD . Line_Form ; 

F0URIER2 . Transform_Type; 

integer: 

integer  :=  -10000; 
integer  : =  1 0000 ; 

float  ;=  0.0; 


float; 

natural; 

character; 

constant  integer  :=  255; 


READ_SSI_IMAGE(  Name,  Char_Count,  Input  ); 

—  [2]  Convert  input  to  complex  image  array 

for  X  in  Input ' first!  1) .. Ir^t' last(  1) 
loop 

for  Y  in  Ir^t'first(2) .  .Input 'last(2) 
loop 

if  Ir¥>ut{X,Y)  =  Bdge_Value  then 
In_Work(X,Y) .Real_N  :=  10.0; 

else 

In_Work(X,Y)  .Real_N  :=  0.0; 
end  if; 

In_Work(X,Y)  .Ijng_N  :=  0.0; 
end  loop; 
end  loop; 

—  [3]  Take  Fourier  Transform 

put_line!"  %REM  -  Taking  2D  Fourier  Transform...  "); 

X_Arg  :=  4.0  *  atan(l.O); 

Y_Arg  ;=  4.0  *  atan(l.O); 

for  X  in  In_Work ' first ( 1 ) . . In_Wbrk ' last ( 1 ) 
loop 

for  Y  in  In_Work'first(2) . .In_Work'last(2) 
loop 

X_Arg_Cos  :=  cos(  float (X)  *  X_Arg  ); 

X_Arg_Sin  .-=  -sin(  float(X)  •  X_Arg  ); 

Y_Arg_Cos  :=  cos(  float (Y)  •  Y_Arg  ); 

Y_Arg_Sin  :=  -sin(  float (Y)  •  Y_Arg  ); 

Real_Shift  ;=  X_Arg_Cos  •  Y_Arg_Cos  -  X_Arg_Sin  *  Y_Arg_Sin; 
Iing_Shift  :=  X_Arg_Sin  *  Y_Arg_Sin  +  X_Arg_Cos  *  Y_Arg_Cos; 

In_Real  :=  In_Work(X,Y) .Real_N; 

In_Img  ;=  In_Work(X,Y) .Img_N; 

In_Work(X,Y) .Real_N  ;=  In_Real  *  Real_Shift  -  In_Img  ♦  ljng_Shift; 
In_Work(X, Y) . Iing_N  :=  In_Img  *  Real_Shift  +  In_Real  *  Img_Shift; 

end  loop; 
end  loop; 

Direction  ;=  forward; 

'IVio_DFT(  In_Work,  Direction)  ; 

X_Atg  ;=  4.0  ♦  atan(l.O); 

Y_Arg  :=  4.0  *  atan(l.O) ; 
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for  X  in  In_Work ' first ( 1 ) . . In_Work ' last ( 1 ) 
loop 

for  Y  in  In_Work ' first ( 2 ) . . In_Work ' last ( 2 ) 
loop 


X_Arg_Cos 

X_Arg_Sin 

Y_Arg_Cos 

Y_Arg_Sin 


=  cos(  float(X)  *  X_Arg  ) 
=  sin(  float (X)  *  X_Arg  ) 
=  cos{  float (Y)  *  Y_Arg  ) 
=  sin(  float (Y)  ♦  Y_Arg  ) 


Real_Shift  :=  X_Arg_Cos  *  Y_Arg__Cos  -  X_Arg_Sin  *  Y_Arg_Sin: 
Img_Shift  :=  X_Arg_Sin  *  Y_Arg__Sin  +  X_Arg_Cos  *  Y_Arg_Cos; 

In_Real  :=  In_Work(X,Y) .Real_N; 

In_Img  : =  In_Work ( X , Y) . Img_N ; 


In_Work(X,Y)  .Real_N 
In_Work(X,Y) .rmg_N 


:=  In_Real  *  Real_Shift  -  In_Img  *  Iing_Shift; 
:=  In_Img  *  Real_Shift  +  In_Real  *  Img_Shift; 


end  loop; 
end  loop; 


new_line; 

put("  SAVE  2D  FOURIER  TRANSFORM?  (y/n)  >  ");  get ( Answer ) ; 

case  Answer  is 
v^ien  'y'  |  'Y'  => 

for  X  in  FT_Output '  first  { 1 ) , .  Fr__Output '  last  ( 1 ) 
loop 

for  Y  in  ET^Output ' first ( 2 ) . . Fr_Output ' last ( 2 ) 
loop 

Temp  :=  float(  In_Work(X,Y) .Real_N**2  +  In_Work(X,Y) . Iing_N**2 
Fr_Output(X,Y)  :=  integer(  sqrt(Temp)  +  0.5  ); 

if  Fr_Output(X,Y)  >  Max_Output  then 
Max_Output  :=  FT_Output(X,Y) ; 
end  if; 

if  5T_0utput(X,Y)  <  Min_Output  then 
Min_Output  : =  Fr_0utput ( X , Y ) ; 
end  if; 

end  loop; 
end  loop; 

for  X  in  FT_0utput ' first ( 1 ) . . FT_0utput ' last ( 1 ) 
loop 

for  Y  in  E*r_0utput ' first ( 2 ) . ,  rr_0utput ' last ( 2 ) 
loop 
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Tanp  ;=  255.0  *  ( float {FT_Output(X,Y )-Min_Output) )  /  float ( Max_Output )  ; 
FT_Output(X,Y)  ;=  integer!  Tanp  ); 
end  loop; 
end  loop: 

Max_Output  : =  - 10000 ; 

Min_0utput  ;=  10000; 

new_line; 

put_line("  ENTER  NAME  OF  FT  FILE" ) ; 

GET_FILENAME(  Name,  Char_Count  ) ; 

SAVE_SSI_IMAGE(  Name,  Char_Count,  F.r_0utput); 

when  others  => 

mill ; 

end  case; 

—  [4]  Convert  to  Polar  Coordinates 

put_line("  SJREM  -  Converting  to  polar  coordinates...”); 

— RECT_T0_P0LAR(  In_Work,  Out_Wbrk  ); 

X_Center  :=  (  In_Work ' last ( 1 )  -  In_Work' first! 1)  )  /  2; 

Y_Center  :=  !  In_Work'last!2)  -  InJWork' first U)  )  /  2; 

for  Theta  in  0ut_Work ' first ! 1 ) . . Out_Work ' last ! 1 ) 
loop 

for  Rho  in  Out_Work'first!2) . ,0ut_Work'last!2) 
loop 

0ut_Work ! Theta, Rho ) .Real_N  ;=  0.0; 

Out_Work!Theta,Pttio) . Img_N  :=  0.0; 

Nearest_X  :=  integer!  float!Rho)  *  cosd! float! Theta ) )  +  0.5  )  +  X_Center: 
Nearest_Y  ;=  integer!  float! Rho)  *  sind! float! Theta) )  +0.5  )  +  Y_Center; 

if  Nearest_X  >=  In_Wbrk ' first ! 1 )  and  Nearest_X  <=  In_Work ' last ! 1 )  and 
Neatrest_Y  >=  In_Work'first!2)  and  Nearest_Y  <=  In_Work'last!2)  then 

0ut_Work ! Theta , Rho ) . Real_N  : =  In_Work! Nearest_X , 257-Nearest_Y ) . Real_N ; 
0ut_Work ! Theta, Rho ). Img_N  ;=  In_Work!Nearest_X, 257-Necu:est_Y) . Img_N; 

end  if; 

end  loop; 
end  loop: 

new_line; 

put!”  SAVE  POLAR  2D  FOURIER  TRANSFORM?  !y/n)  >  ");  get ! Answer ) ; 
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case  Answer  is 
when  'y'  |  'Y'  => 


for  Theta  in  FT_RT_Output ' f irst ( 1 ) . .  Fr_RT_Output ' last ( 1 ) 
loop 

for  Rho  in  FT_RT_Output ' first { 2 ) . .  Fr_RT_Output ' last { 2 ) 
loop 

Temp  ;=  float(  Out_Work(Theta,Rho) .Real_N**2 

+  cXit_Work(Theta,Rho) . Img_N**2  ); 
rr_Rr_Output( Theta, Rho)  :=  integer(  sqrt(Temp)  +  0.5  ) 

if  Fr_RT_0utput( Theta, Rho)  >  Max_Output  then 
Max_Output  :=  Fr_RT_0utput( Theta, Rho) ; 
end  if; 

if  FT_RT_0utput{ Theta, Rho)  <  Min_Output  then 
Min_Output  :=  FT_RT_0utput( Theta, Rho) ; 
end  if; 

end  loop; 
end  loop; 

for  Theta  in  FT_RT_0utput ' first (1 ) . . FT_RT_Output ' last ( 1 ) 
loop 

for  Rho  in  FT_RT_CXitput ' first ( 2 ) . .  ET_RT_0utput ' last ( 2 ) 
loop 

Tenp  :=  255.0  •  {  float(  BT_RT_Output ( Theta , Rho )  -  Min 
/  f  loat  ( Max_0utput ) ; 

FT_RT_0utput( Theta, Rho)  :=  integer (  Temp  ); 
end  loop; 
end  loop; 

Max_0utput  :=  -10000; 

Min_Output  :=  10000; 

new_line; 

put_line("  ENTER  NAME  OF  POLAR  FT  FILE"); 

GET_FILENAME(  Name,  Char_Count  ) ; 

SAVE_SSI_IMAGE(  Name,  Char_Count,  Fr_RT_Output ) ; 

when  others  => 

nxill ; 

end  case; 

—  [5]  Take  inverse  Fourier  Transform  along  the  Rho  direction 

put_line("  S5REM  -  Taking  radial  inverse  Fourier  Transform..."); 

for  Rho  in  Rho_Work ' first . . Rho_Work ' last 
loop 


Rho_Work(Rho) .ReaI_N  ;=  0.0: 

Rho_Work(Rho) . Img_N  :=  0.0; 
end  lcx>p; 

Direction  :=  inverse; 

for  Theta  in  0ut_Work ' first ( 1 ) . . Out_Work ' last ( 1 ) 
loop 

—  for  Rho  in  Out_Work ' first ( 2 ) . . Out_Work ' last ( 2 ) 
for  Rho  in  0. . 127 

loop 

Rho_Work  ( Rho ) .  Real_N  :  =  OutJWork  ( Theta ,  Rho+1 ) .  Real_N ; 

Rho_Wbrk(Ey:io)  . Iing_N  ;=  Out_Work(Theta,Rho+l)  . Img_N; 
end  loop; 

One_DFT(  Rho_Work,  Direction  ) ; 

for  Rho  in  Output 'first( 2) . .Output 'last(2) 

—  for  Rho  in  0..127 
loop 

Tenip_0ut  Rho_Wbrk(255-Rho)  .Real_N**2  +  Rho_Work(255-Rho)  .Iaig_N**2; 
Ten?3_0ut  :=  Rho_Wbrk(Rho) .Real_N**2  +  Rho_Work(Rho).Img_N*»2; 

Output (Theta, Rho)  ;=  integer{  sqrt { Ten¥)_0ut )  +  0.5  ); 

if  Max_Output  <  Output (Theta, Rho)  then 
Max_0utput  :=  Output (Theta, Rho) ; 
end  if; 

if  Min_Output  >  Output  (Theta, Rho)  then 
Min_0utput  :*  Output (Theta, Rho) ; 
end  if; 

end  loop; 

for  Rho  in  129. ,181 

—  loop 

—  Output  (Theta,  Rho)  :=  0: 
end  loop; 

end  loop; 

new_line; 

putC'Max  Output  =  '');  put(Max_Output,3,10) ;  new_line; 
put ("Min  Output  =  ");  put ( Min_Output , 3 , 10 ) ; 

—  [6]  Save  in  2D  Radon  file 

for  Theta  in  Output' first( 1) . .Output 'last(l) 
loop 

for  Rho  in  Output ’ first(2) . .Output ' last(2) 
loop 

Temp  ;=  float (  Output ( Theta , Rho )  -  Min  Output  )  *  255.0; 


Temp  :=  Temp  /  float  {  Max_Output  -  Min__Output  ); 

Out^t( Theta, Rho)  :=  integer(  Tonp  ); 
end  loop; 
end  loop; 

—  [6.1]  form  histogram 
new_line; 

put("  FORM  HISTOGRAM?  (y/n)  >  ");  get  (Answer);  new_line: 

case  Answer  is 
when  'y'  |  'Y'  => 

HIST0_GRAM1_2D(  Output,  Output_Histogram  ); 

for  Index  in  Output_Histogram' first . .Output_Histogram' last 
loop 

Output_Histogram_f loat ( Index )  :=  float (  Output_Histogram ( Index )  ) 
end  loop; 

new_line: 

put_line("  ENTER  NAME  OF  MATRIXX  HISTOGRAM  FILE"); 

GET_FILENAME(  Name,  Char_Count  ); 

MATRIX_X_E»L0T1D(  Name,  Char_Count,  Output_Histogram_float  ); 
when  others  =>  null; 
end  case; 


new_line(2 ) ; 

put_line("  ENTER  RADON  SSI  FILENAME"); 
GET_FILENAME(  Name,  Char_Count  ) ; 
SAVE_SSI_IMAGE(  Name,  Char_Ccfunt,  Output  ) 

new_line(2) : 

put_line("  END  FOURIER_RADON  ••*"); 
end  FOURIER  RADON; 
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INV  FR 


—  Purpose;  To  generate  the  orginal  image  from  the  Radon  Transform  throvigh 

—  a  radial  ID  Fourier  Transform  followed  by  a  polar  to  reactangulair  coordinate 

—  transformation  and  an  inverse  2D  Fourier  Transform. 


with  DATA_STANDARD: 
with  DATA_CONVERISONS: 
with  F0URIER2; 
with  SSI_IO; 
with  FILE_IO; 
with  text_io; 
with  integer_text_io ; 
with  float  math  lib; 


use  DATA_STANDARD; 
use  DATA_CONVERSIONS; 
use  P0URIER2; 
use  SSI_IO; 
use  FILE_IO; 
use  text_io; 
iise  integer_text_io ; 
use  float  math  lib; 


procedure  INV_FR  is 

Input,  RFT_Output 

Output,  Fr_Output 

In_Work 

Out_WQrk 

Rho_Work 

Output_Histogram 

Output_Histogram_float 

Name 

Direction 

Nearest_Rho, 

Nearest_Theta 

Block_DC 

X_Center,  Y_Center 

Max_Output 

Min_Output 

Temp 

Temp_Out , 

pi 

X_Arg_Cos,  X_Arg_Sin. 

Y_Arg_Cos,  Y_Arg_Sin, 

In_Real,  In_Img, 

Real_Shift,  Img_Shift 

XI,  Y1 

Char_Cciunt 

Answer 

Bdge_Valtie 

begin  —  INV_FR 


DATA_STANDARD . Image_Array_2d (  0 . . 359 ,  0 . . 181  ) ; 
DATA_STANPARD . Image_Array_2d (  1 . . 256 ,  1 , . 256  ) ; 
DATA_STANnARD.Con55lex_Array_R2d(  0..359,  0..127  ) 
DATA_STANDARD. Complex  Array_R2d(  1..256,  1..256  ) 
DATA  STANDARD. Ccn5ile!xlArray_Rld(  0..255  ); 
DATA2sTANDARD.Image_Array_ld(  0..255  ); 

DATA  STANDARD. Array_ld{  0..255  ) ; 

DATa'sTANDARD . Line_Form ; 

F0URIER2 . Transf orm_Type ; 


integer: 

constant  integer 
constant  integer 


integer 

integer 

float 


=  -10000; 
=  10000; 
=  0.0; 


constant  float  :=  4.0  •  atan(l.O): 


float ; 
float; 
natural ; 
character; 

constant  integer  ;=  255; 


—  [1]  Get  input  image 
new_line: 

put_line(''  ENTER  RADON  INPUT  IMAGE"); 

GET_FILENAME(  Name,  Char_Caunt  ) : 

READ_SSI_IMAGE (  Name ,  Char_Count ,  Input  ) ; 

—  [2]  Convert  input  to  ccn^lex  image  array  and  take  Radial  Fourier  Transform 
put_line("  5SREM  -  Taking  radial  Fourier  Transform..."); 

Direction  :=  forward; 

for  Theta  in  InJWork ' f irst { 1 ) . . In_Vk3rk ' last ( 1 ) 
loop 

for  Rho  in  In_Work ' first ( 2 ) . . InJWork ' last ( 2 ) 
loop 

Y_Arg_Cos  :=  cos(  float(Rho)  *  pi  ) ; 

Y_Arg_Sin  ;=  -sin(  float (Rho)  *  pi  ) ; 

Y_Arg_Cos  ;=  1.0; 

Y~Arg_SIn  :=  1.0; 

Rho_Work(Rho) .Real_N  :=  float (Input (Theta, Rho ) )  *  Y_Arg_Cos; 

Rho_Work ( Rho ) . Img_N  : *  float ( Input { Theta , Rho ) )  *  Y_Arg_Sin ; 

end  loop; 

for  Rho  in  In_Work'last(2)+l. .Rho_Wbrk'last 
loop 

Rho_Work(Rho) .Real_N  :=  0.0; 

Rho_Work(Rho) .Img_N  :=  0.0; 
end  loop; 

0ne_DFr(  Rho_Work,  Direction  ) ; 

for  Rho  in  In_Work'first(2) . .In_Work'last(2) 
loop 

In_Vfork( Theta, Rho) .Real_N  ;=  Rho_Work(Rho) .Real_N; 

In_Work ( Theta , Eho ) . Img_N  : =  Rho_Work ( Rho ) . Img_N ; 
end  loop; 

end  loop; 

new_line; 

put("  SAVE  RADIAL  FOURIER  TRANSFORM?  (y/n)  >  ");  get ( Answer ) ; 


case  Answer  is 
when  'y'  1  'Y'  => 


for  Theta  in  RFT_Output ' first ( 1 ) . . RFT_Oatput ' last ( 1 ) 
loop 

for  Rho  in  RFT_Output ' first ( 2 ) . . In_Work ' last ( 2 ) 
loop 

Temp  : =  float {  In_Work ( Theta , Rho ) , Real_N*  *  2 
+  In_Work ( Theta , Rho ) . Img_N*  *2  ) ; 

RFT_Output( Theta, Rho)  :=  integer(  sqrt(Temp)  +  0.5  ) 

if  RFr_Output( Theta, Rho)  >  Max_Output  then 
Max_Output  :=  RFT_Output (Theta, Rho ) ; 
end  if; 

if  RFT_Output( Theta, Rho)  <  Min_Output  then 
Min_Outpit  ;=  RFT_Output{ Theta, Rho) ; 
end  if; 

end  loop; 

for  Rho  in  In_Work ' last ( 2 )  +1  . .  RFT_Output ' last ( 2 ) 
loop 

RFr_Output( Theta, Rho)  :=  Min_Output; 
end  loop; 

end  loop; 

for  Theta  in  RET_Output ' first { 1 ) . . RFT_Output ' last ( 1 ) 
loop 

for  Rho  in  RFT_Output '  first  ( 2 ) . . RBT_Output '  last  ( 2 ) 
loop 

Temp  :=  255.0  *  (  float(  RFT_Output( Theta, Rho) 

-  Min_0utput  )  )  /  float (Max_0utput ) ; 
F?FT_Output( Theta, Rho)  :=  integer (  Ten?)  ); 
end  loop; 
end  loop; 

Max_Output  :=  -10000; 

Min_Output  :=  10000; 

new_line; 

put_line("  ENTER  NAME  OF  RFT  FILE"); 

GEr_FILENAME(  Nanve,  Char_Count  ); 

SAVE_SSI_IMAGE(  Name,  Char_Count,  RFT_0utput); 

when  others  => 

mill ; 


end  case; 

—  [3]  Polar  to  Cartesian 


for  X  in  Out_Wbrk'first(l) . .CXit_Work'last(l) 

IcxDp 

for  Y  in  OutJWbrk ' f irst ( 2 ) . . Out_Work ' last ( 2 ) 
loop 

XI  :=  float (X  -  X_Center); 

Y1  :=  float (257  -  Y  -  Y_Center) ; 

Nearest_Rho  :=  integer(  sqrt{  Xl**2  +  Yl**2  )  ) ; 

if  XI  =  0.0  and  Y1  >=  0.0  then 
Nearest_Theta  :=  90; 
elsif  XI  =  0.0  and  VI  <  0.0  then 
Nearest_Theta  :=  -90; 
else 

Nearest_Theta  :=  integer (  atan2d(  Yl,  XI  )  ); 
end  if; 

if  Nearest_Theta  <  0  then 

Nearest_Theta  :=  Nearest_Theta  +  360; 
end  if; 

if  Nearest_Theta  in  In_Work'range(l)  aM  Nearest_Rho  in  In_Work ' range ( 2 ) 
then 

0ut_Work(X,Y) .Real_N  :=  In_Work(Nearest_Theta,Nearest_Rho).Real_N; 
Out_Work(X,Y) .Img_N  :=  InJ*3rk(Nearest_Theta,Nearest_Rho) .Inig_N; 
end  if; 

end  loop; 
end  loop; 

new_line; 

put("  SAVE  RECT  FOURIER  TRANSFORM?  (y/n)  >  ");  get { Answer ) ; 

case  Answer  is 
when  'y'  |  'Y'  => 

for  X  in  Fr_0utput ' first ( 1 ) . .  ET_0utput ' last ( 1 ) 
loop 

for  Y  in  FT_0utput ' first ( 2 ) . .  ET^Output ' last { 2 ) 
loop 

Tanp  :=  float(  Gut_Work{X,Y) .Real_N**2  +  0ut_W3rk(X, Y) . Img_N**2  ); 
Fr_0utput(X,Y)  :=  integer(  sqrt(Temp)  +  0.5  ) ; 

if  FT_Output(X,Y)  >  Max_Output  then 
Max_0utput  ;=  BT_0utput(X,Y) ; 
end  if; 

if  FT_0utput(X,Y)  <  Min_Output  then 
Min_0utput  ; =  FT_0utput ( X , Y ) ; 
end  if; 
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end  loop; 
end  loop: 

for  X  in  ET_Output ' first ( 1 ) . .  FT^Output ’ last ( 1 ) 
loop 

for  Y  in  Fr_Output '  first  ( 2 ) .  .  FT_Output '  last  { 2 ) 
loop 

Ten^)  :=  255.0  *  (  float(  ET_Output(X,Y)  -  Min_0utput  )  ) 
/  float ( Max_0utput ) ; 

Fr_0utput(X,Y)  :=  integer(  Temp  ); 
end  loop; 
end  loop; 

Max_0utput  :=  -10000; 

Min_Output  :=  10000; 

new_line; 

put_line("  ENTER  NAME  OF  FT  FILE"); 

GET_FILENAME(  Name,  Char_Count  ); 

SAVE_SSI_IMAGE(  Name,  C3Tar_Count,  Fr_0utput); 

when  others  => 

null ; 

end  case; 

—  [4]  Take  2D  inverse  Fourier  Transform 

put_line("  ajREM  -  Taking  2D  Fourier  Transform...  "); 

for  X  in  Out_Work ' first ( 1 ) . . 0ut_W3rk ' last { 1 ) 
loop 

for  Y  in  Out_Work ' f irst ( 2 ) . . OutJWork ' last ( 2 ) 
loop 


X_Arg_Cos 

X_Arg_Sin 

Y_Arg_Cos 

Y_Arg_Sin 


=  cos(  float(X)  *  pi  ) ; 

=  -sin(  float(X)  *  pi  ); 
=  cos(  float(Y)  ♦  pi  ) ; 

=  -sin(  float (Y)  *  pi  ) ; 


Real_Shift  :=  X_Arg__Cos  *  Y_Arg_Cos  -  X_Arg_Sin  *  Y_Arg_Sin; 
Img_Shift  :=  X_Arg_Sin  *  Y_Arg_Sin  +  X_Arg_Cos  *  Y_Arg_Cos; 

In_Real  ;=  Out_Work(X,Y) .Real_N; 

In_Img  ; =  0ut_Work ( X , Y ) , Img_N ; 

0ut_Work(X,Y) .Real_N  :=  In_Real  •  Real_Shift  -  In_Iing  *  Img_Shift; 
0ut_Work(X,Y)  .Img_N  :=  In_Iing  *  Real_Shift  +  In_Real  *  Img__Shift; 


end  loop; 
end  loop; 


Direction  :=  inverse; 

Two_DET(  Out_Work,  Direction) ; 

—  [4.1]  Block  DC  component 

for  X  in  X_Center  -  Block_DC  . .  X_Center  +  Block_DC 

^°°for  Y  in  Y_Center  -  Block_DC  . .  Y_Center  +  Block_DC 
loop 

Out_Work{X,Y) .Real_N  :=0.0: 

Out_Work{X,Y)  .Iing_N  :=  0.0; 
end  loop; 
end  loop; 

for  X  in  Output ' f iret ( 1 ) . . Output ' last ( 1 } 
loop 

for  Y  in  Output ' first ( 2 ) . . Output ' last ( 2 ) 
loop 

Temp  :=  float(  Out_Wbrk(X,Y) .Real_N**2  +  OutJfork(X,Y) .  Img_N**2  ); 
Output(X,Y)  :=  integer(  sqrt{Tanp)  +  0.5  ); 

if  Output (X,Y)  >  Max_Output  then 
Max_0utput  :=  Output(X,Y) ; 
end  if; 

if  0utput{X,Y)  <  Min_0utput.  then 
Min_0utput  :=  Output (X,Y) ; 
end  if; 


end  loop; 
end  loop; 

for  X  in  Output ' first ( 1 ) . . Output ' last ( 1 ) 

loop  ,  ,  , 

for  Y  in  Output ' first ( 2 ) . . Output ' last { 2 ) 

^°°Temp  :=  255.0  *  (  float (  Output (X,Y)  -  Min_0utput  )  )  /  float (Max_0utput) 
Output (X,Y)  :=  integer (  Tanp  ); 
end  loop; 
end  loop; 

—  [5]  form  histogram 

new  line;  ,  ,  .  .  ^ 

put("  FOIW  HISTOGRAM?  (y/n)  >  ");  get (Answer);  new_line; 

case  Answer  is 
when  'y'  |  'Y'  => 

HISTO  GHAMl  2D(  Output,  Output_Histogram  ); 


for  Index  in  Output_Histogram' first . .Output_Histogram' last 
loop 

CXitput_Histogram_fIoat( Index)  :=  float (  Output_Histogram{ Index) 
end  loop; 

new_line; 

put_line("  ENTER  NAME  OF  MATRIXX  HISTOGRAM  FILE"); 

GET_FILENAME(  Name,  Chsur_Count  ) ; 

MATRIX_X_PL0T1D(  Name,  Char_Count,  Output_Histogram_float  ); 
when  others  =>  null; 
end  case; 

—  [6]  Save  image 
new_line( 2) ; 

put_line("  ENTER  SSI  FILENAME"); 

GET_FILENAME(  Name,  Char_Caunt  ); 

SAVE_SSI_IMAGE(  Name,  Char_Count,  Output  ); 

new_line(2) ; 

put_line("  END  INV_FR  ***"); 
end  INV  FR; 


RAD0N2 


Purpose:  Perform  2D  Radon  Transform  on  input  image 


With  DATA_STANDARD; 
with  FILE_IO; 
with  SSI_IO; 
with  float_math_lib: 
with  text  io; 


use  DATA_STANDARD; 
use  FILE_IO; 
use  SSI_IO; 
vise  float_math_lib; 
use  text  io; 


with  integer_text_io ;  use  integer_text_io; 


DATA__STANDARD .  Image_Array_2d  (  1 . .  256 ,  1 . .  256  ) ; 

DATA_STANDARD . Image_Array_2d (  0..359,  0..181  ); 

float ; 

integer ; 

integer  :=  0; 

integer  ;=  2000; 

integer  :=  255; 

integer  :=  128; 

DATA_STANDARD . Line_Form ; 
natural; 


procedure  RAD0N2  is 

Ir^t 

Output 

Theta_Float,  X.  Y 

Rho_0ut,  Tenp_0ut 

Max_0utput 

Min_0utput 

Edge_Value 

Shift 

Name 

Char_Count 

begin  —  RAD0N2 

—  [  1  ]  Get  ir^TUt  image 


new_line; 

put_line("  NAME  OF  SSI  IMAGE  ?"); 
GET_FILENAME(Name,  Char_Count); 
READ_SSI_IMAGE(Name,  Char_Count,  Input); 

—  [2]  Set  Radon  space  to  zero 
new_line; 

put_line("  3JF?EM  -  Ccanputing  Radon  Transform..."); 

for  Theta  in  Output ' first ( 1 ). .Output ' last( 1 ) 
loop 

for  Rho  in  Output ' first ( 2) . .Output 'last( 2) 
loop 

Output(Theta,  Rho)  :=  0; 
end  loop; 
end  loop; 

—  [3]  Take  Radon  Transform 


new  line: 


for  XI  in  Input '  first  ( 1 ) . .  Ir^t '  last  ( 1 ) 
loop 

for  Y1  in  Input ' first ( 2) . .Input 'last{ 2) 
loop 

if  Input (XI, Yl)  >  0  then 

for  Theta  in  Output ' first { 1 Output ' last ( 1 ) 
loop 

X  :=  float {XI  -  Shift): 

Y  ;=  float(257  -  Yl  -  Shift); 

Theta_Float  :=  float ( Theta ) ; 

Rho_Out  :=  integer (  X  *  coi^(Theta_Float)  +  Y  *  sind ( Theta_Float )  ) 
if  Rho_Out  >  0  then 

Output  ( Theta,  Rho_Out)  :=  Output  (Theta,  Rho_Out)  +  I:^t(Xl,Yl); 
elsif  Rho_Out  =  0  and  Theta  <  180  then 

Output  ('Iheta,Rho_Out)  :=  Output  (Theta, Rho_0ut)  +  Input  ( XI ,  Yl )  ; 
end  if; 

end  loop: 
end  if; 
end  loop; 
end  loop; 

new_line: 

for  Theta  in  Output 'first(l ). .Output 'last (1) 
loop 

for  Rho  in  Output ' first (2) . .Output ' last (2) 
loop 

if  Max_0utput  <  Output (Theta, Rho)  then 
Max_Output  :=  Output (Theta, Rho) ; 
end  if; 

if  Min_0utput  >  Output (Theta, Rho)  then 
Min_0utput  : =  Output ( Theta , Rho ) ; 
end  if; 

end  loop; 
end  loop; 

new_line; 

for  Theta  in  Output ' f irst( 1) . .Output 'last(l) 
loop 

for  Rho  in  Output ' first { 2) . .Output 'last( 2) 
loop 

Output  ( Theta ,  Rho ) 
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THRESHOLD  HS  IMAGE  (SSI  -  Format) 


Purpose;  Threshold  Hoxjgh  Transforro  Image 


with  DATA_STANDARD; 
with  FILE_IO: 
with  SSI_IO; 
with  text  io; 


use  DATA_STANDARD; 
vise  FILE_IO: 
use  SSI_IO; 
use  text  io; 


with  integer_text_io ;  use  integer_text_io ; 


procedure  THRESHOLD_HS  is 

Input 

Name 

Lcwer,  Upper,  Max 

Char_Count 

Answer 

begin  —  THRESHOLD_HS 
—  [1]  Get  image 


new_line; 

put_line("  NAME  OF  THE  INPUT  SSI  FILE?"); 
GET_FILENAME(NsBne,  Char_Count) ; 
READ_SSI_IMAGE  ( Name ,  Char_Count ,  Input ) ; 

for  Theta  in  Ir^t'first(l) .  .Input 'last(l) 
loop 

for  Rho  in  Input 'first( 2) . .Input 'last(2) 
loop 

if  Max  <  Ir^t(Theta,  Rho)  then 
Max  ;=  Input(Theta,  Rho) ; 
end  if; 
end  loop; 
end  loop: 

new_line; 

put ("MAXIMUM  VALUE  =  ");  put(Max) ;  new_line: 

—  [2]  Get  threshold  value 

begin 

GET_LCWER: 

loop 

new  line; 


DATA_STANDARD . Image_ArraY_2d ( 0 . . 359 , 0 . , 1 81 ) ; 

DATA_STANDARD . Line_Form ; 

integer  range  0 . . 255 ; 

natural; 

character : 


Jigs' 


.*  'y~. 
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si  IsIvXlv 


put_line("  ENTER  LOWER  BOUND  (0  -  255}") 
get(Lciwer) ; 

put("  LOWER  BOUND  =  ") ; 
put { Lcwer ) ; 
put("?  (Y/N)  =>  "); 
get (Answer) ; 


case  Answer  is 

when  'Y'  |  'y'  => 
exit  GET_LCWER; 
when  others  => 
null; 
end  case; 


end  loop  GET_L0WER; 


GET_UPPER: 

loop 

new_line; 

put_line("  ENTER  UPPER  BOUND  (0  -  255)"); 
get ( Upper ) ; 

put("  UPPER  BOUND  =  "); 
put (Upper) ; 
put("?  (Y/N)  =>  "); 
get ( Answer ) ; 


case  Answer  is 


exit  GET_UPPER; 
when  others  ®> 
null; 
end  case; 


end  loop  GET_UPPER; 


exception 

when  Data_Error  => 

put_line( "  Invalid  entry. " ) ; 

put_line("  Enter  positive  integer  valxoe  (0  -  255).") 

end; 


—  [3]  Threshold  and  save  in  SSI  format 

put_line( "Begin  thresholding. . . " ) ; 

for  Rho  in  Input ' first ( 2) . .Input 'last(2) 
loop 

for  Theta  in  Input ' f irst(l) . .Input 'last a) 
loop 

if  Input (Theta, Rho)  >  Upper  or  Input ( Theta , Rho )  <  Lower  then 
Input (Theta, Rho)  :=  0; 
end  if; 

•  end  loop ; 
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end  loop; 
new_line: 

put_line(''  NAME  OF  OUTPUT  SSI  FILE?"); 
GET_FILENAME(Name,  Char_Covint )  ; 
SAVE_SSI_IMAGE{Nanie,  Char_Count,  Input), 

new_line; 

put_line("  ***  END  THFESHOLD_HS  ***"); 
new__line: 

end  THRESHOLD  HS; 


_ ^m**************^************************************************************* 


CXJCLUDE 


—  Purpose:  To  block  a  section  of  input  image 


With  DATA_STANDARD; 

with  FILE_IO; 

with  SSI_IO; 

with  text_io: 

with  integer_text_io ; 

with  float  math  lib; 


use  DATA_STANDARD; 

use  FILE_IO; 

use  SSI_IO; 

use  text_io; 

use  integer_text_io; 

lise  float  math  lib; 


procedure  OCCLUDE  is 

Input 

Name 

Char_Count 

Lower_Left_X,  Lawer_Left_Y 
Box_Size 

X_Center,  Y_Center 

X_Float,  Y_Float 

Radius 

Answer 

Choice 


nATA_STANDARD . Image_Array_2d (  1 . . 256 ,  1 . . 256  ) ; 

DATA_STANDARD . Line_Form ; 

natural; 

integer ; 

integer ; 

integer ; 

float; 

integer; 

character ; 

integer  r’ange  0 . .  3  ; 


begin  —  OCCLUDE 
—  [1]  Get  ir^t  image 
new_line; 

put_line("  ENTER  NAME  OF  SSI  INPUT"); 
GET_FILENAME{  Name,  Char_Count  ) ; 
READ_SSI_IMAGE(  Name,  Char_Count,  Input  ); 


—  [2]  Get  blocking  parameters 


GET_LOOP: 

loop 

new_line(3)  ; 

put_line("  ***  OCCLUSION  TYPES  **•"); 
new_line; 

put("  [1]  Circular");  new_line(2); 

put("  [2]  Square");  new_line(2); 

put("  ENTER  CHOICE  >  ");  get (Choice);  new_line; 

case  Choice  is 

when  1  I  2  =>  exit  GET_LOOP; 


when  others  =>  null; 
end  case; 

end  loop  GET_LOOP; 

case  Choice  is 
v«rfien  1  => 

CIRC: 

loop 

put_line{"  ***  Circle  Parameters  ***"): 

put("  [1]  Center  =  (");  put(X_Center,4, 10) ; 

put('’,");  put(Y_Center,4, 10) ;  put(")");  new_line(2); 

put("  [2]  Circle  Radius  =  ");  put (Radius, 3, 10) ;  new_line(2); 

put("  [0]  Continue");  new_line(2); 

put("  Enter  Choice  =>  ");  get(Choice);  new_line; 

case  Choice  is 
when  1  => 

CENTER_X; 

loop 

new_line: 

put_line("  ENTER  X  Value  (-127.. 128)  *>  ");  get(X_Center) ; 

put("  X  value  =  ");  put{X_Center.4, 10) ;  put("  (Y/N)?");  get{Answer) 

case  Answer  is 

when  'y'  |  'Y'  =>  exit  CENTER_X; 
when  others  =>  null; 
end  case; 

end  loop  CENTER_X; 

CENTER_Y: 

loop 

new_line; 

put_line("  ENTER  Y  Value  (-127.. 128)  =>  ");  get { Y_Center ) ; 

put{"  Y  value  =  ");  put { Y_Center , 4 , 10 ) ;  put("  (Y/N)?");  get{Answer) 

case  Answer  is 

when  'y'  |  'Y'  =>  exit  CENTER_Y; 
v«*)en  others  =>  null; 
end  case; 

end  loop  CENTER_Y; 

when  2  => 

G1ET_RADIUS: 

loop 


new_lane; 

put_line(''  Enter  Radi^ls  Value  (1,.128)  =>  ");  get ( Radius); 
put("  Radius  =  ");  put(Radi\is,3, 10) ;  put("  (Y/N)?");  get(Answer) 

case  Answer  is 

'y'  I  'Y'  =>  exit  GET_RADIUS; 
when  others  =>  null; 
end  case; 

end  loop  GET_RADIUS; 
when  others  =>  null; 
end  case; 
end  loop  CIRC; 

put_line("  %REM  -  Setting  valioes  to  zero..."); 

for  X  in  Input' first (1) . .Input' last (1) 
loop 

for  Y  in  Input ' first { 2) . .Input 'last (2) 
loop 

X_Float  ;=  (  float(X)  -  128.0  )  -  float ( X_Center ) ; 

Y_Float  ;=  {  257.0  -  float(Y)  -  128.0  )  -  float(Y_Center) ;  . 

if  sqrt(  X_Float**2  +  Y_Float*»2  )  <=  float (Radius)  then 
Input (X,Y)  :=  0; 
end  if; 

end  loop; 
end  loop; 

when  2  => 


put_line("  Box  Parameters  ***"); 

put("  [1]  Lower  Left  Comer  =  {");  put(Lower_Left_X,4,10); 
put (",");  put(Lcwer_Left_Y,4,10) ;  put(")"); 
put(")");  new_line(2); 

put("  [2]  Size  =  ");  put(Box_Size,3, 10) ;  new_line(2); 

put("  [0]  Continue");  new_line(2); 

put("  Enter  Choice  =>  ");  get{Choice);  new_line; 

case  Choice  is 
when  1  => 

LCMER_LEFT: 

loop 


new_line; 

put_line("  Enter  X  value  of  lower  left  comer  (-127..  128) 
get(Lcwer_Left_X)  ; 

put_line("  Enter  Y  valioe  lower  left  comer  (-127..  128)  => 
get(Lower_Left_Y) ; 

put("  Y  valiie  =  ");  pait(Lcwer_Left_Y,4, 10) ;  new_line; 
put("  X  value  =  ");  put(Lower_Left_X,4, 10) ;  put("  (Y/N)?") 

get(Answer);  new_line; 

case  Answer  is 

when  'y'  |  'Y'  =>  exit  LOWER_LEFT; 
when  others  =>  null; 
end  case; 

end  loop  LCWER_LEFT; 

when  2  => 

SIZE; 

loop 

new_line; 

put_line("  Enter  size  of  box  (1..255)  =>  ");  get(Box_Size) 
put("  Box  size  =  ");  put(Box_Size,3, 10) ;  put("  (Y/N)?"); 
get(Answer) ; 

case  Answer  is 
when  'y'  |  'Y'  =>  exit  SIZE; 
vAien  others  =>  null; 
end  case; 

end  loop  SIZE; 

when  others  => 

exit  BOX; 

end  case; 

end  loop  BOX; 

put_line("  ajREM  -  Setting  values  to  zero..."); 

for  X  in  Input ' first( 1) . .Input 'last(l) 
loop 

for  Y  in  Input '  first  ( 2 ) . .  Ir^xit '  last  { 2 ) 
loop 

if  X  -  128  >  Lower  Left  X 
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and  X  -  128  <  Lower_Left_X  +  Box_Size  then 

if  256  -  Y  -  128  >  Lcwer_Left_Y 

and  256  -  Y  -  128  <  Lower_Left_Y  +  Box_Size  then 

Ir^JUt(X,Y)  ;=  0; 

end  if; 

end  If; 

end  loop; 
end  loop; 

v*ien  others  =>  null; 
end  case; 
new_line(2)  ; 

put_line("  ENTER  NAME  OF  SSI  EDGE  FILE")  ; 

GET_FILENAME(  Name,  Char_Count  ) ; 

SAVE_SSI_IMAGE(  Name,  Char_Count,  Input  ); 

new_line(2) ; 

put_line("  ***  END  OCCLUDE  ***"): 
new_line(2); 


end  OCCLUDE 


_ m*ritmti*m*****«**m*************m************************************************ 

— 

ROTATE_EDGE  IMAGE  (SSI  -  Format) 

—  Purpose:  Rotate  an 

image  by  theta 

with  DATA  STANDARD; 

use  DATA  STANDARD; 

witii  FILE_IO; 

use  FILE  10; 

with  SSI_IO; 

use  SSI  10; 

with  te?rt_io; 

xise  text  io; 

with  integer,. text  io; 

use  integer_text_io; 

with  fIoatjnath_lib; 

vise  float_inath_llb; 

procedure  ROTATE_EDGE  is 

Input,  Output 

:  DATA  STANDARD. Image  Array  2d(l. .256,1. .256) ; 

Name 

;  DATA_STANDARD.Line_Form; 

Bc3ge_VaIue 

:  constant  integer  :=  255; 

X_Center,  Y_Center 

:  constant  float  :=  128.0; 

Angle 

:  integer  range  0 . . 359 ; 

X.  Y,  Radiias, 

Theta  In,  Theta  Out 

;  float; 

X_Out,  Y_Out 

:  integer; 

Char  Count 

:  natural; 

Answer 

:  chau^cter; 

begin  —  ROTATE_EDGE 

—  [1]  Get  image 

new  line; 

put  line("  NAME  OF  THE 

INPUT  SSI  FILE?”); 

GET  FILENAME  ((feune,  Char 

Count) ; 

READ_SSI_IMAGE{Naine,  Char_Count,  Ii^nit); 

—  [2]  Set  output  image 

to  zero 

for  Y1  in  Output ' first ( 2 ) . . Output ' last ( 2 ) 

loop 

for  XI  in  Output ' first ( 1 ) . . Output ' last ( 1 ) 

loop 

Output (XI, Yl)  :=  ' 

0; 

end  loop; 

end  loop; 

—  [2]  Get  rotation  angle 

X 


GET_ANC3LE: 

loop 

new_line; 

putC  ENTER  ROTATION  ANGLE  (0  ..  359)  >  ");  get(Angle);  new_line; 
put("  Angle  =  ");  ^t(Angle,  3,  10);  put("?  (Y/N)  =>  "); 
get  (Answer);  new_line; 

case  Answer  Is 
when  'Y'  \  'y'  => 
exit  GET_ANGLE; 
when  others  => 
null; 
end  case; 

end  loop  GET_ANGLE; 

—  [3]  Shift  and  save  in  SSI  format 

put_line("  %REM  -  Rotating  edge  image..."); 

for  Y1  in  Input'first(2) .  .Ir^t'last(2) 
loop 

for  XI  in  Input ' first ( 1 ) . . Input ' last ( 1 ) 
loop 

if  Ir53ut{Xl,yi)  =  Bdge_Value  then 

X  :»  floatvXl)  -  X  Center; 

Y  :=  float (Yl)  -  y'Center; 

Radius  ;»  sqrt(X**2  +  Y**2); 

Theta_In  :=  atan2d(y,X) ; 

Theta_0ut  :*  'Iheta_In  +  float  (Angle ) ; 

X_0ut  :=  integer  (Radius  •  cosd(Theta_Out)  +  X_Center) ; 

Y_0ut  ;*«  integer  (Radius  •  sind(Theta_0ut)  +  y_Center)  ; 

if  X_0ut  in  Output' ranged)  and  Y_0ut  in  Output 'range( 2)  then 

Output ( X_0ut ,  Y_0ut)  :*  Edge_Value; 

end  if; 

end  if; 
end  loop; 
end  loop; 

new_line; 

put_line("  NAME  OF  OUTPUT  SSI  FILE?"): 
jET_FrLENAI«(NamB,  Char_Count); 

SAVE_SSI_IMAGE(N«ae,  Char_Count,  Output); 


MOVE_EDGE  IMAGE  (SSI  -  Format) 


—  Purpose:  Translate  an  ecJge  image  by  input  x  and  y  increments 


With  DATA_STANDARD; 
with  FILE_IO; 
with  SSI_IO; 
with  text_io; 
with  integer_text_io : 


use  nATA_STANDARD; 
use  FILE_IO; 
use  SSI_IO; 
use  text_io; 

\ase  integer_text_io; 


procedure  MOVE_EDGE  is 


Input,  Output 
N^e 

X_Shift,  Y_Shift 
Edge_Value 
Chetr_Count 
Answer 

begin  —  MOVE_EDGE 


nATA_STANnARD . Image_Array_2d ( 1 . . 256 , 1 . . 256 ) ; 

DATA_STANDARD . Line_Form ; 

integer  range  -128.. 128; 

constant  integer  :=  255; 

natux^I; 

character ; 


—  [1]  Get  image 
new_liiie; 

put_line("  NAME  OF  THE  INPUT  SSI  FILE?"); 
GET_FILQlAME(Name,  Chau:_Count ) ; 
READ_SSI_IMAGE(Nanie,  Char_Count,  Ir^t); 

—  [2]  Set  output  image  to  zero 

for  Y  in  Output ' first { 2) . .Output 'last (2) 
loop 

for  X  in  Output 'first(l) . .Output 'last(l) 
loop 

Output(X,Y)  :=  0; 
end  loop; 
end  loop; 


—  [2]  Get  shift  values 


begin 

GET_X_SHIFT: 

loop 

new_line; 

put_line("  ENTER  SHIFT  IN  X  (-128  ..  128)"); 
get(X_Shift) ; 


put{"  X  SHIFT  =  ") ; 
put(X_Shift) ; 
put("?  (Y/N)  =>  "); 
get(Answer) ; 

case  Answer  is 

whei  'Y'  I  'y'  => 
exit  GET_X_SHIFr; 

Mhen  others  => 
null; 
end  case; 

end  loop  GET_X_SHIFr: 

3ET_Y_SHlFr: 

loop 

nev/_line: 

put_line(''  ENTER  SHIFT  IN  Y  (-128  ..  128)"); 

get(Y_Shift) ; 

put("  Y  SHIFT  =  "); 

put(Y_Shift) ; 

put("?  (Y/N)  =>  "): 

get(Answer) ; 

case  Answer  is 

when  'Y'  1  'y'  => 
exit  GET_Y_SHIFT; 
when  others  => 
null; 
end  case; 

end  loop  GET_Y_SHIFr; 

exception 

v*en  Data_Error  => 

put_line( "  Invalid  entry, " ) ; 

put_line("  Biter  integer  Vetlue  (-128  ..  128)."); 

end; 

—  [3]  Shift  and  save  in  SSI  format 

put_line("  Jtf?EM  -  Shifting  edge  points..."); 

for  Y  in  Input '  first  ( 2 ) . .  Input '  last  ( 2 ) 
loop 

for  X  in  Input ' first( 1) . .Input 'last (1) 
loop 

if  Input  (X,Y)  =  Edge_Value  then 

Output (X  +  X_Shift,  Y  +  Y_Shift)  :=  Edge_Value: 
end  if; 


^Hi***»^**»*******:tt*m:ii*^iii*jitmm***************************************l‘********** 


EXPAND 


—  Purpose;  Perform  linear  e}^)ansion  of  image  for  display 


With  DATA_STANDARD; 

with  FILE_IO: 

with  SSI_IO; 

with  text_io; 

with  integer_text_io; 

with  floatjnath_lib: 

procedure  EXPAND  is 


use  DATA_STANDARD; 

use  FILE_IO; 

use  SSI_IO; 

use  text_io; 

use  integer_text_io; 

use  float  math  lib; 


Input 

Ma^c 

Min 

Name 

Char  Count 


DATA_STANDARD . Image_Array_2d( 0 . . 359 , 0 . ■ 181 ) : 
integer  :=  0; 
integer  :=  300; 

DATA_STANDARD . Line_Form ; 
integer; 


begin  —  EXPAND 
new_line; 

put_line( "NAME  OF  HS  SSI  INPOT  FILE?"); 

GET  FILENAME (Name,  Char  Count) ; 

READ_SSI_IMAGE(IfeBiie,  Char_Count,  Ii^jut); 

new_line; 

put_line(''  SREM  -  Perfcming  linear  e^qjansion, . . ")  ; 

for  XI  in  Input' first( 1) . .Input 'last(l) 
loop 

for  Y1  in  Input 'first (2) . .Input 'last(2) 
loop 

if  Max  <  Input (XI, Yl)  then 
Max  :=  Input(Xl,Yl)  ; 
end  if; 

if  Min  >  Ir^t(Xl,Yl)  then 
Min  :=  Il^t(Xl,Yl)  ; 
end  if; 
end  loop; 
end  loop; 

for  XI  in  Input ' first ( 1 ) . . Input ' last ( 1 ) 
loop 

for  Yl  in  Inpwt '  f iret  ( 2 ) . .  Ir^t '  last  ( 2 ) 
loop 

Input(Xl,Yl)  :*  (  255  *  (  Input(Xl,Yl)  -  Min  )  )  /  (  Max  -  Min  ) ; 


m 


■  ^ 


m 
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na^_line; 

put_line("NAME  OF  SSI  OUTPUT  FILE"); 
GET_FILENAME(Name,  Char_Count) ; 
SAVE_SSI_IMAGE(Name,  Char_Count,  Input) 

new_line; 

put_line("  END  EXPAND  ***"); 
nfiw_line; 
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